From 321ee0b1bc863730f2c52e281654a0b73eecfc01 Mon Sep 17 00:00:00 2001 From: Marshall Ward Date: Tue, 25 Oct 2022 13:49:13 -0400 Subject: [PATCH] Generic I/O layer for infra and netCDF This patch introduces a generalized I/O interface, through a new class, `MOM_file`, for reading and writing of axis-based model fields. A common API for interacting with files is defined in the `MOM_file` class, and two implementations are provided: * `MOM_infra_file`, using the "infra" framework (i.e. FMS) * `MOM_netcdf_file`, using native netCDF framework This separation allows us to define certain generic functions by class, and platform-specific functions by type. It will also allow for removal of legacy FMS1 operations which are no longer used for the majority of I/O but are still required for isolated MOM files. This change will allow us to move incompatible files from the FMS to the netCDF implementation, while still preserving a common structure for both files. The majority of the details of `MOM_file` and its subclasses are defined in `MOM_io_file.F90`, which is exclusively accessed through `MOM_io.F90`. The netCDF implementation, designed to be used by `MOM_netcdf_file` but is designed as a standalone system, is defined in `MOM_netcdf.F90`. *Interface* `MOM_file` includes the following functions: * `open` * `close` * `flush` * `register_axis` * `register_field` * `write_attribute` * `write_field` * `file_is_open` * `get_file_info` * `get_file_fields` * `get_field_atts` * `read_field_chksum` Most are designed to resemble the existing MOM I/O operations. Note that some of these have not yet been implemented for `MOM_netcdf_file`, since they are never used in the model. `MOM_file_infra` includes the following additional operations: * `get_file_times` * `get_file_fieldtypes` See documentation for usage of these functions. The "axis"/"field" model from FMS1 has been preserved, where axes are associated with variables which contain the grid point, and fields are defined with respect to the file's axes. Operations such as field counters will exclude any variables associated with the axes. The `axistype` and `fieldtype` from FMS have been replaced with new abstractions (`MOM_axis` and `MOM_field`). Internally, these point to either the FMS types or equivalent netCDF types. *Implementation* Each file type contains an instance of its native type, as well as lists of associated axes and fields. List are implemented as linked lists of names (stored as `label`) which are used to identify and extract or write the axis/field to the internal type. The mechanics of this are largely hidden from the user and can be changed in the future, if needed, without disruption to the rest of the codebase. The current netCDF implementation very closely mirrors the FMS infra behavior, but this can be relaxed or modified as needed. netCDF I/O is currently only designed for serial I/O, with all of the data on the root PE. Further development would be needed to support any kind of parallel I/O. *Current Usage* Two functions have been transferred from infra to the netCDF I/O: * `Depth_list.nc` * `ocean.stats.nc` The following legacy functions and types from FMS have been preserved: * `create_file` * `reopen_file` * `file_type` * `open_file` * `get_file_info` * `get_file_fields` * `get_file_times` This is primarily to preserve compatibility with SIS2, but may also be useful for other code, such as model drivers. These may be phased out in the future, however. --- config_src/infra/FMS2/MOM_io_infra.F90 | 4 +- src/ALE/MOM_hybgen_regrid.F90 | 13 +- src/ALE/MOM_regridding.F90 | 14 +- src/diagnostics/MOM_sum_output.F90 | 72 +- src/framework/MOM_io.F90 | 442 ++++- src/framework/MOM_io_file.F90 | 1654 +++++++++++++++++ src/framework/MOM_netcdf.F90 | 755 ++++++++ src/framework/MOM_restart.F90 | 48 +- src/ice_shelf/MOM_ice_shelf.F90 | 2 +- .../MOM_coord_initialization.F90 | 12 +- .../MOM_shared_initialization.F90 | 12 +- 11 files changed, 2854 insertions(+), 174 deletions(-) create mode 100644 src/framework/MOM_io_file.F90 create mode 100644 src/framework/MOM_netcdf.F90 diff --git a/config_src/infra/FMS2/MOM_io_infra.F90 b/config_src/infra/FMS2/MOM_io_infra.F90 index c49c124ae0..54b9dfb78b 100644 --- a/config_src/infra/FMS2/MOM_io_infra.F90 +++ b/config_src/infra/FMS2/MOM_io_infra.F90 @@ -509,8 +509,8 @@ subroutine open_ASCII_file(unit, file, action, threading, fileset) ! This checks if open() failed but did not raise a runtime error. inquire(unit, opened=is_open) if (.not. is_open) & - call MOM_error(FATAL, 'open_ASCII_file: File ' // trim(filename) // & - ' failed to open.') + call MOM_error(FATAL, & + 'open_ASCII_file: File "' // trim(filename) // '" failed to open.') ! NOTE: There are two possible mpp_write_meta functions in FMS1: ! - call mpp_write_meta( unit, 'filename', cval=mpp_file(unit)%name) diff --git a/src/ALE/MOM_hybgen_regrid.F90 b/src/ALE/MOM_hybgen_regrid.F90 index 271caa7ad6..f89e15d930 100644 --- a/src/ALE/MOM_hybgen_regrid.F90 +++ b/src/ALE/MOM_hybgen_regrid.F90 @@ -7,7 +7,8 @@ module MOM_hybgen_regrid use MOM_EOS, only : EOS_type, calculate_density use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, assert use MOM_file_parser, only : get_param, param_file_type, log_param -use MOM_io, only : close_file, create_file, file_type, fieldtype, file_exists +use MOM_io, only : create_MOM_file, file_exists +use MOM_io, only : MOM_infra_file, MOM_field use MOM_io, only : MOM_read_data, MOM_write_field, vardesc, var_desc, SINGLE_FILE use MOM_string_functions, only : slasher use MOM_unit_scaling, only : unit_scale_type @@ -210,20 +211,20 @@ subroutine write_Hybgen_coord_file(GV, CS, filepath) character(len=*), intent(in) :: filepath !< The full path to the file to write ! Local variables type(vardesc) :: vars(3) - type(fieldtype) :: fields(3) - type(file_type) :: IO_handle ! The I/O handle of the fileset + type(MOM_field) :: fields(3) + type(MOM_infra_file) :: IO_handle ! The I/O handle of the fileset vars(1) = var_desc("dp0", "meter", "Deep z-level minimum thicknesses for Hybgen", '1', 'L', '1') vars(2) = var_desc("ds0", "meter", "Shallow z-level minimum thicknesses for Hybgen", '1', 'L', '1') vars(3) = var_desc("Rho_tgt", "kg m-3", "Target coordinate potential densities for Hybgen", '1', 'L', '1') - call create_file(IO_handle, trim(filepath), vars, 3, fields, SINGLE_FILE, GV=GV) + call create_MOM_file(IO_handle, trim(filepath), vars, 3, fields, & + SINGLE_FILE, GV=GV) call MOM_write_field(IO_handle, fields(1), CS%dp0k, scale=CS%coord_scale) call MOM_write_field(IO_handle, fields(2), CS%ds0k, scale=CS%coord_scale) call MOM_write_field(IO_handle, fields(3), CS%target_density, scale=CS%Rho_coord_scale) - call close_file(IO_handle) - + call IO_handle%close() end subroutine write_Hybgen_coord_file !> This subroutine deallocates memory in the control structure for the hybgen module diff --git a/src/ALE/MOM_regridding.F90 b/src/ALE/MOM_regridding.F90 index e28f2c5e82..53072909a5 100644 --- a/src/ALE/MOM_regridding.F90 +++ b/src/ALE/MOM_regridding.F90 @@ -6,8 +6,9 @@ module MOM_regridding use MOM_error_handler, only : MOM_error, FATAL, WARNING, assert use MOM_file_parser, only : param_file_type, get_param, log_param use MOM_io, only : file_exists, field_exists, field_size, MOM_read_data -use MOM_io, only : vardesc, var_desc, fieldtype, SINGLE_FILE -use MOM_io, only : create_file, MOM_write_field, close_file, file_type +use MOM_io, only : vardesc, var_desc, SINGLE_FILE +use MOM_io, only : MOM_infra_file, MOM_field +use MOM_io, only : create_MOM_file, MOM_write_field use MOM_io, only : verify_variable_units, slasher use MOM_unit_scaling, only : unit_scale_type use MOM_variables, only : ocean_grid_type, thermo_var_ptrs @@ -2212,8 +2213,8 @@ subroutine write_regrid_file( CS, GV, filepath ) character(len=*), intent(in) :: filepath !< The full path to the file to write type(vardesc) :: vars(2) - type(fieldtype) :: fields(2) - type(file_type) :: IO_handle ! The I/O handle of the fileset + type(MOM_field) :: fields(2) + type(MOM_infra_file) :: IO_handle ! The I/O handle of the fileset real :: ds(GV%ke), dsi(GV%ke+1) if (CS%regridding_scheme == REGRIDDING_HYBGEN) then @@ -2231,10 +2232,11 @@ subroutine write_regrid_file( CS, GV, filepath ) vars(2) = var_desc('ds_interface', getCoordinateUnits( CS ), & 'Layer Center Coordinate Separation', '1', 'i', '1') - call create_file(IO_handle, trim(filepath), vars, 2, fields, SINGLE_FILE, GV=GV) + call create_MOM_file(IO_handle, trim(filepath), vars, 2, fields, & + SINGLE_FILE, GV=GV) call MOM_write_field(IO_handle, fields(1), ds) call MOM_write_field(IO_handle, fields(2), dsi) - call close_file(IO_handle) + call IO_handle%close() end subroutine write_regrid_file diff --git a/src/diagnostics/MOM_sum_output.F90 b/src/diagnostics/MOM_sum_output.F90 index 7797a266dd..6dbe4997d6 100644 --- a/src/diagnostics/MOM_sum_output.F90 +++ b/src/diagnostics/MOM_sum_output.F90 @@ -13,8 +13,9 @@ module MOM_sum_output use MOM_forcing_type, only : forcing use MOM_grid, only : ocean_grid_type use MOM_interface_heights, only : find_eta -use MOM_io, only : create_file, file_type, fieldtype, flush_file, reopen_file, close_file -use MOM_io, only : file_exists, slasher, vardesc, var_desc, write_field, MOM_write_field +use MOM_io, only : create_MOM_file, reopen_MOM_file +use MOM_io, only : MOM_infra_file, MOM_netcdf_file, MOM_field +use MOM_io, only : file_exists, slasher, vardesc, var_desc, MOM_write_field use MOM_io, only : field_size, read_variable, read_attribute, open_ASCII_file, stdout use MOM_io, only : axis_info, set_axis_info, delete_axis_info, get_filename_appendix use MOM_io, only : attribute_info, set_attribute_info, delete_attribute_info @@ -125,9 +126,9 @@ module MOM_sum_output !! to stdout when the energy files are written. integer :: previous_calls = 0 !< The number of times write_energy has been called. integer :: prev_n = 0 !< The value of n from the last call. - type(file_type) :: fileenergy_nc !< The file handle for the netCDF version of the energy file. + type(MOM_netcdf_file) :: fileenergy_nc !< The file handle for the netCDF version of the energy file. integer :: fileenergy_ascii !< The unit number of the ascii version of the energy file. - type(fieldtype), dimension(NUM_FIELDS+MAX_FIELDS_) :: & + type(MOM_field), dimension(NUM_FIELDS+MAX_FIELDS_) :: & fields !< fieldtype variables for the output fields. character(len=200) :: energyfile !< The name of the energy file with path. end type sum_output_CS @@ -603,13 +604,11 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci energypath_nc = trim(CS%energyfile) // ".nc" if (day > CS%Start_time) then - call reopen_file(CS%fileenergy_nc, trim(energypath_nc), vars, & - num_nc_fields, CS%fields, SINGLE_FILE, CS%timeunit, & - G=G, GV=GV) + call reopen_MOM_file(CS%fileenergy_nc, trim(energypath_nc), vars, & + num_nc_fields, CS%fields, SINGLE_FILE, CS%timeunit, G=G, GV=GV) else - call create_file(CS%fileenergy_nc, trim(energypath_nc), vars, & - num_nc_fields, CS%fields, SINGLE_FILE, CS%timeunit, & - G=G, GV=GV) + call create_MOM_file(CS%fileenergy_nc, trim(energypath_nc), vars, & + num_nc_fields, CS%fields, SINGLE_FILE, CS%timeunit, G=G, GV=GV) endif endif @@ -863,35 +862,35 @@ subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, dt_forci endif endif - call write_field(CS%fileenergy_nc, CS%fields(1), real(CS%ntrunc), reday) - call write_field(CS%fileenergy_nc, CS%fields(2), toten, reday) - call write_field(CS%fileenergy_nc, CS%fields(3), PE, reday) - call write_field(CS%fileenergy_nc, CS%fields(4), KE, reday) - call write_field(CS%fileenergy_nc, CS%fields(5), H_0APE, reday) - call write_field(CS%fileenergy_nc, CS%fields(6), mass_lay, reday) - - call write_field(CS%fileenergy_nc, CS%fields(7), mass_tot, reday) - call write_field(CS%fileenergy_nc, CS%fields(8), mass_chg, reday) - call write_field(CS%fileenergy_nc, CS%fields(9), mass_anom, reday) - call write_field(CS%fileenergy_nc, CS%fields(10), max_CFL(1), reday) - call write_field(CS%fileenergy_nc, CS%fields(11), max_CFL(2), reday) + call CS%fileenergy_nc%write_field(CS%fields(1), real(CS%ntrunc), reday) + call CS%fileenergy_nc%write_field(CS%fields(2), toten, reday) + call CS%fileenergy_nc%write_field(CS%fields(3), PE, reday) + call CS%fileenergy_nc%write_field(CS%fields(4), KE, reday) + call CS%fileenergy_nc%write_field(CS%fields(5), H_0APE, reday) + call CS%fileenergy_nc%write_field(CS%fields(6), mass_lay, reday) + + call CS%fileenergy_nc%write_field(CS%fields(7), mass_tot, reday) + call CS%fileenergy_nc%write_field(CS%fields(8), mass_chg, reday) + call CS%fileenergy_nc%write_field(CS%fields(9), mass_anom, reday) + call CS%fileenergy_nc%write_field(CS%fields(10), max_CFL(1), reday) + call CS%fileenergy_nc%write_field(CS%fields(11), max_CFL(2), reday) if (CS%use_temperature) then - call write_field(CS%fileenergy_nc, CS%fields(12), 0.001*Salt, reday) - call write_field(CS%fileenergy_nc, CS%fields(13), 0.001*salt_chg, reday) - call write_field(CS%fileenergy_nc, CS%fields(14), 0.001*salt_anom, reday) - call write_field(CS%fileenergy_nc, CS%fields(15), Heat, reday) - call write_field(CS%fileenergy_nc, CS%fields(16), heat_chg, reday) - call write_field(CS%fileenergy_nc, CS%fields(17), heat_anom, reday) + call CS%fileenergy_nc%write_field(CS%fields(12), 0.001*Salt, reday) + call CS%fileenergy_nc%write_field(CS%fields(13), 0.001*salt_chg, reday) + call CS%fileenergy_nc%write_field(CS%fields(14), 0.001*salt_anom, reday) + call CS%fileenergy_nc%write_field(CS%fields(15), Heat, reday) + call CS%fileenergy_nc%write_field(CS%fields(16), heat_chg, reday) + call CS%fileenergy_nc%write_field(CS%fields(17), heat_anom, reday) do m=1,nTr_stocks - call write_field(CS%fileenergy_nc, CS%fields(17+m), Tr_stocks(m), reday) + call CS%fileenergy_nc%write_field(CS%fields(17+m), Tr_stocks(m), reday) enddo else do m=1,nTr_stocks - call write_field(CS%fileenergy_nc, CS%fields(11+m), Tr_stocks(m), reday) + call CS%fileenergy_nc%write_field(CS%fields(11+m), Tr_stocks(m), reday) enddo endif - call flush_file(CS%fileenergy_nc) + call CS%fileenergy_nc%flush() if (is_NaN(En_mass)) then call MOM_error(FATAL, "write_energy : NaNs in total model energy forced model termination.") @@ -1233,13 +1232,13 @@ subroutine write_depth_list(G, US, DL, filename) ! Local variables type(vardesc), dimension(:), allocatable :: & vars ! Types that described the staggering and metadata for the fields - type(fieldtype), dimension(:), allocatable :: & + type(MOM_field), dimension(:), allocatable :: & fields ! Types with metadata about the variables that will be written type(axis_info), dimension(:), allocatable :: & extra_axes ! Descriptors for extra axes that might be used type(attribute_info), dimension(:), allocatable :: & global_atts ! Global attributes and their values - type(file_type) :: IO_handle ! The I/O handle of the fileset + type(MOM_netcdf_file) :: IO_handle ! The I/O handle of the fileset character(len=16) :: depth_chksum, area_chksum ! All ranks are required to compute the global checksum @@ -1259,8 +1258,8 @@ subroutine write_depth_list(G, US, DL, filename) call set_attribute_info(global_atts(1), depth_chksum_attr, depth_chksum) call set_attribute_info(global_atts(2), area_chksum_attr, area_chksum) - call create_file(IO_handle, filename, vars, 3, fields, SINGLE_FILE, extra_axes=extra_axes, & - global_atts=global_atts) + call create_MOM_file(IO_handle, filename, vars, 3, fields, SINGLE_FILE, & + extra_axes=extra_axes, global_atts=global_atts) call MOM_write_field(IO_handle, fields(1), DL%depth, scale=US%Z_to_m) call MOM_write_field(IO_handle, fields(2), DL%area, scale=US%L_to_m**2) call MOM_write_field(IO_handle, fields(3), DL%vol_below, scale=US%Z_to_m*US%L_to_m**2) @@ -1268,8 +1267,7 @@ subroutine write_depth_list(G, US, DL, filename) call delete_axis_info(extra_axes) call delete_attribute_info(global_atts) deallocate(vars, extra_axes, fields, global_atts) - call close_file(IO_handle) - + call IO_handle%close() end subroutine write_depth_list !> This subroutine reads in the depth list from the specified file diff --git a/src/framework/MOM_io.F90 b/src/framework/MOM_io.F90 index e8df89b268..6c2dc9df34 100644 --- a/src/framework/MOM_io.F90 +++ b/src/framework/MOM_io.F90 @@ -15,15 +15,17 @@ module MOM_io use MOM_io_infra, only : read_field, read_vector use MOM_io_infra, only : read_data => read_field ! Deprecated use MOM_io_infra, only : read_field_chksum -use MOM_io_infra, only : file_type, file_exists, get_file_info, get_file_fields -use MOM_io_infra, only : open_file, open_ASCII_file, close_file, flush_file, file_is_open -use MOM_io_infra, only : get_field_size, fieldtype, field_exists, get_field_atts -use MOM_io_infra, only : get_file_times, axistype, get_axis_data, get_filename_suffix -use MOM_io_infra, only : write_field, write_metadata, write_version +use MOM_io_infra, only : file_exists +use MOM_io_infra, only : open_ASCII_file, close_file, file_is_open +use MOM_io_infra, only : get_field_size, field_exists, get_field_atts +use MOM_io_infra, only : get_axis_data, get_filename_suffix +use MOM_io_infra, only : write_version use MOM_io_infra, only : MOM_namelist_file, check_namelist_error, io_infra_init, io_infra_end use MOM_io_infra, only : APPEND_FILE, ASCII_FILE, MULTIPLE, NETCDF_FILE, OVERWRITE_FILE use MOM_io_infra, only : READONLY_FILE, SINGLE_FILE, WRITEONLY_FILE use MOM_io_infra, only : CENTER, CORNER, NORTH_FACE, EAST_FACE +use MOM_io_file, only : MOM_file, MOM_infra_file, MOM_netcdf_file +use MOM_io_file, only : MOM_axis, MOM_field use MOM_string_functions, only : lowercase, slasher use MOM_verticalGrid, only : verticalGrid_type @@ -33,21 +35,33 @@ module MOM_io use netcdf, only : NF90_strerror, NF90_inquire_dimension use netcdf, only : NF90_NOWRITE, NF90_NOERR, NF90_GLOBAL, NF90_ENOTATT, NF90_CHAR +! The following are not used in MOM6, but may be used by externals (e.g. SIS2). +use MOM_io_infra, only : axistype ! still used but soon to be nuked +use MOM_io_infra, only : fieldtype +use MOM_io_infra, only : file_type +use MOM_io_infra, only : get_file_info +use MOM_io_infra, only : get_file_fields +use MOM_io_infra, only : get_file_times +use MOM_io_infra, only : open_file +use MOM_io_infra, only : write_field + implicit none ; private ! These interfaces are actually implemented in this file. -public :: create_file, reopen_file, cmor_long_std, ensembler, MOM_io_init +public :: create_MOM_file, reopen_MOM_file, cmor_long_std, ensembler, MOM_io_init +public :: MOM_field public :: MOM_write_field, var_desc, modify_vardesc, query_vardesc, position_from_horgrid public :: open_namelist_file, check_namelist_error, check_nml_error public :: get_var_sizes, verify_variable_units, num_timelevels, read_variable, read_attribute public :: open_file_to_read, close_file_to_read ! The following are simple pass throughs of routines from MOM_io_infra or other modules. -public :: file_exists, open_file, open_ASCII_file, close_file, flush_file, file_type -public :: get_file_info, field_exists, get_file_fields, get_file_times, get_filename_appendix +public :: file_exists, open_ASCII_file, close_file +public :: MOM_file, MOM_infra_file, MOM_netcdf_file +public :: field_exists, get_filename_appendix public :: fieldtype, field_size, get_field_atts public :: axistype, get_axis_data public :: MOM_read_data, MOM_read_vector, read_field_chksum -public :: slasher, write_field, write_version_number +public :: slasher, write_version_number public :: io_infra_init, io_infra_end public :: stdout_if_root public :: get_var_axes_info @@ -67,6 +81,15 @@ module MOM_io !> These encoding constants are used to indicate the discretization position of a variable public :: CENTER, CORNER, NORTH_FACE, EAST_FACE +! The following are not used in MOM6, but may be used by externals (e.g. SIS2). +public :: create_file +public :: reopen_file +public :: file_type +public :: open_file +public :: get_file_info +public :: get_file_fields +public :: get_file_times + !> Read a field from file using the infrastructure I/O. interface MOM_read_data module procedure MOM_read_data_0d @@ -87,6 +110,11 @@ module MOM_io !> Write a registered field to an output file, potentially with rotation interface MOM_write_field + module procedure MOM_write_field_legacy_4d + module procedure MOM_write_field_legacy_3d + module procedure MOM_write_field_legacy_2d + module procedure MOM_write_field_legacy_1d + module procedure MOM_write_field_legacy_0d module procedure MOM_write_field_4d module procedure MOM_write_field_3d module procedure MOM_write_field_2d @@ -147,23 +175,70 @@ module MOM_io character(len=:), allocatable :: att_val !< The values of this attribute end type attribute_info - integer, public :: stdout = stdout_iso !< standard output unit integer, public :: stderr = stderr_iso !< standard output unit contains -!> Routine creates a new NetCDF file. It also sets up fieldtype -!! structures that describe this file and variables that will -!! later be written to this file. -subroutine create_file(IO_handle, filename, vars, novars, fields, threading, timeunit, & - G, dG, GV, checksums, extra_axes, global_atts) - type(file_type), intent(inout) :: IO_handle !< Handle for a files or fileset that is to be +!> `create_MOM_file` wrapper for the legacy file handle, `file_type`. +!! NOTE: This function may be removed in a future release. +subroutine create_file(IO_handle, filename, vars, novars, fields, threading, & + timeunit, G, dG, GV, checksums, extra_axes, global_atts) + type(file_type), intent(inout) :: IO_handle + !< Handle for a files or fileset that is to be opened or reopened for + !! writing + character(len=*), intent(in) :: filename + !< full path to the file to create + type(vardesc), intent(in) :: vars(:) + !< structures describing fields written to filename + integer, intent(in) :: novars + !< number of fields written to filename + type(fieldtype), intent(inout) :: fields(:) + !< array of fieldtypes for each variable + integer, optional, intent(in) :: threading + !< SINGLE_FILE or MULTIPLE + real, optional, intent(in) :: timeunit + !< length of the units for time [s]. The default value is 86400.0, for 1 + !! day. + 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. + integer(kind=int64), optional, intent(in) :: checksums(:,:) + !< checksums of vars + type(axis_info), optional, intent(in) :: extra_axes(:) + !< Types with information about some axes that might be used in this file + type(attribute_info), optional, intent(in) :: global_atts(:) + !< Global attributes to write to this file + + type(MOM_infra_file) :: new_file + type(MOM_field) :: new_fields(novars) + + new_file%handle_infra = IO_handle + + call create_MOM_file(new_file, filename, vars, novars, new_fields, & + threading=threading, timeunit=timeunit, G=G, dG=dG, GV=GV, & + checksums=checksums, extra_axes=extra_axes, global_atts=global_atts) + + IO_handle = new_file%handle_infra + call new_file%get_file_fieldtypes(fields(:novars)) +end subroutine create_file + + +!! Create a new netCDF file and register the MOM_fields to be written. +subroutine create_MOM_file(IO_handle, filename, vars, novars, fields, & + threading, timeunit, G, dG, GV, checksums, extra_axes, global_atts) + class(MOM_file), intent(inout) :: IO_handle !< Handle for a files or fileset that is to be !! opened or reopened for writing character(len=*), intent(in) :: filename !< full path to the file to create type(vardesc), intent(in) :: vars(:) !< structures describing fields written to filename integer, intent(in) :: novars !< number of fields written to filename - type(fieldtype), intent(inout) :: fields(:) !< array of fieldtypes for each variable + type(MOM_field), intent(inout) :: fields(:) !< array of fieldtypes for each variable integer, optional, intent(in) :: threading !< SINGLE_FILE or MULTIPLE real, optional, intent(in) :: timeunit !< length of the units for time [s]. The !! default value is 86400.0, for 1 day. @@ -186,10 +261,10 @@ subroutine create_file(IO_handle, filename, vars, novars, fields, threading, tim logical :: use_layer, use_int, use_periodic logical :: one_file, domain_set, dim_found logical, dimension(:), allocatable :: use_extra_axis - type(axistype) :: axis_lath, axis_latq, axis_lonh, axis_lonq - type(axistype) :: axis_layer, axis_int, axis_time, axis_periodic - type(axistype), dimension(:), allocatable :: more_axes ! Axes generated from extra_axes - type(axistype) :: axes(5) ! The axes of a variable + type(MOM_axis) :: axis_lath, axis_latq, axis_lonh, axis_lonq + type(MOM_axis) :: axis_layer, axis_int, axis_time, axis_periodic + type(MOM_axis), dimension(:), allocatable :: more_axes ! Axes generated from extra_axes + type(MOM_axis) :: axes(5) ! The axes of a variable type(MOM_domain_type), pointer :: Domain => NULL() type(domain1d) :: x_domain, y_domain integer :: position, numaxes, pack, thread, k, n, m @@ -244,9 +319,9 @@ subroutine create_file(IO_handle, filename, vars, novars, fields, threading, tim if (domain_set) one_file = (thread == SINGLE_FILE) if (one_file) then - call open_file(IO_handle, filename, OVERWRITE_FILE, threading=thread) + call IO_handle%open(filename, action=OVERWRITE_FILE, threading=thread) else - call open_file(IO_handle, filename, OVERWRITE_FILE, MOM_domain=Domain) + call IO_handle%open(filename, action=OVERWRITE_FILE, MOM_domain=Domain) endif ! Define the coordinates. @@ -326,28 +401,23 @@ subroutine create_file(IO_handle, filename, vars, novars, fields, threading, tim "create_file: A vertical grid type is required to create a file with a vertical coordinate.") if (use_lath) & - call write_metadata(IO_handle, axis_lath, name="lath", units=y_axis_units, longname="Latitude", & + axis_lath = IO_handle%register_axis("lath", units=y_axis_units, longname="Latitude", & cartesian='Y', domain=y_domain, data=gridLatT(jsg:jeg)) - if (use_lonh) & - call write_metadata(IO_handle, axis_lonh, name="lonh", units=x_axis_units, longname="Longitude", & + axis_lonh = IO_handle%register_axis("lonh", units=x_axis_units, longname="Longitude", & cartesian='X', domain=x_domain, data=gridLonT(isg:ieg)) - if (use_latq) & - call write_metadata(IO_handle, axis_latq, name="latq", units=y_axis_units, longname="Latitude", & + axis_latq = IO_handle%register_axis("latq", units=y_axis_units, longname="Latitude", & 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", & + axis_lonq = IO_handle%register_axis("lonq", units=x_axis_units, longname="Longitude", & 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), & + axis_layer = IO_handle%register_axis("Layer", units=trim(GV%zAxisUnits), & longname="Layer "//trim(GV%zAxisLongName), cartesian='Z', & sense=1, data=GV%sLayer(1:GV%ke)) - if (use_int) & - call write_metadata(IO_handle, axis_int, name="Interface", units=trim(GV%zAxisUnits), & + axis_int = IO_handle%register_axis("Interface", units=trim(GV%zAxisUnits), & longname="Interface "//trim(GV%zAxisLongName), cartesian='Z', & sense=1, data=GV%sInterface(1:GV%ke+1)) @@ -367,9 +437,9 @@ subroutine create_file(IO_handle, filename, vars, novars, fields, threading, tim write(time_units,'(es8.2," s")') timeunit endif - call write_metadata(IO_handle, axis_time, name="Time", units=time_units, longname="Time", cartesian='T') + axis_time = IO_handle%register_axis("Time", units=time_units, longname="Time", cartesian='T') else - call write_metadata(IO_handle, axis_time, name="Time", units="days", longname="Time", cartesian='T') + axis_time = IO_handle%register_axis("Time", units="days", longname="Time", cartesian='T') endif ; endif if (use_periodic) then @@ -378,24 +448,24 @@ subroutine create_file(IO_handle, filename, vars, novars, fields, threading, tim ! Define a periodic axis with unit labels. allocate(axis_val(num_periods)) do k=1,num_periods ; axis_val(k) = real(k) ; enddo - call write_metadata(IO_handle, axis_periodic, name="Period", units="nondimensional", & - longname="Periods for cyclical variables", cartesian='T', data=axis_val) + axis_periodic = IO_handle%register_axis("Period", units="nondimensional", & + longname="Periods for cyclical variables", cartesian='T', data=axis_val) deallocate(axis_val) endif do m=1,num_extra_dims ; if (use_extra_axis(m)) then if (allocated(extra_axes(m)%ax_data)) then - call write_metadata(IO_handle, more_axes(m), name=extra_axes(m)%name, units=extra_axes(m)%units, & + more_axes(m) = IO_handle%register_axis(extra_axes(m)%name, units=extra_axes(m)%units, & longname=extra_axes(m)%longname, cartesian=extra_axes(m)%cartesian, & sense=extra_axes(m)%sense, data=extra_axes(m)%ax_data) elseif (trim(extra_axes(m)%cartesian) == "T") then - call write_metadata(IO_handle, more_axes(m), name=extra_axes(m)%name, units=extra_axes(m)%units, & + more_axes(m) = IO_handle%register_axis(extra_axes(m)%name, units=extra_axes(m)%units, & longname=extra_axes(m)%longname, cartesian=extra_axes(m)%cartesian) else ! FMS requires that non-time axes have variables that label their values, even if they are trivial. allocate (axis_val(extra_axes(m)%ax_size)) do k=1,extra_axes(m)%ax_size ; axis_val(k) = real(k) ; enddo - call write_metadata(IO_handle, more_axes(m), name=extra_axes(m)%name, units=extra_axes(m)%units, & + more_axes(m) = IO_handle%register_axis(extra_axes(m)%name, units=extra_axes(m)%units, & longname=extra_axes(m)%longname, cartesian=extra_axes(m)%cartesian, & sense=extra_axes(m)%sense, data=axis_val) deallocate(axis_val) @@ -457,10 +527,10 @@ subroutine create_file(IO_handle, filename, vars, novars, fields, threading, tim pack = 1 if (present(checksums)) then - call write_metadata(IO_handle, fields(k), axes(1:numaxes), vars(k)%name, vars(k)%units, & + fields(k) = IO_handle%register_field(axes(1:numaxes), vars(k)%name, vars(k)%units, & vars(k)%longname, pack=pack, checksum=checksums(k,:)) else - call write_metadata(IO_handle, fields(k), axes(1:numaxes), vars(k)%name, vars(k)%units, & + fields(k) = IO_handle%register_field(axes(1:numaxes), vars(k)%name, vars(k)%units, & vars(k)%longname, pack=pack) endif enddo @@ -468,41 +538,92 @@ subroutine create_file(IO_handle, filename, vars, novars, fields, threading, tim if (present(global_atts)) then do n=1,size(global_atts) if (allocated(global_atts(n)%name) .and. allocated(global_atts(n)%att_val)) & - call write_metadata(IO_handle, global_atts(n)%name, global_atts(n)%att_val) + call IO_handle%write_attribute(global_atts(n)%name, global_atts(n)%att_val) enddo endif - ! Now actualy write the variables with the axis label values - if (use_lath) call write_field(IO_handle, axis_lath) - if (use_latq) call write_field(IO_handle, axis_latq) - if (use_lonh) call write_field(IO_handle, axis_lonh) - if (use_lonq) call write_field(IO_handle, axis_lonq) - if (use_layer) call write_field(IO_handle, axis_layer) - if (use_int) call write_field(IO_handle, axis_int) - if (use_periodic) call write_field(IO_handle, axis_periodic) + ! Now write the variables with the axis label values + if (use_lath) call IO_handle%write_field(axis_lath) + if (use_latq) call IO_handle%write_field(axis_latq) + if (use_lonh) call IO_handle%write_field(axis_lonh) + if (use_lonq) call IO_handle%write_field(axis_lonq) + if (use_layer) call IO_handle%write_field(axis_layer) + if (use_int) call IO_handle%write_field(axis_int) + if (use_periodic) call IO_handle%write_field(axis_periodic) do m=1,num_extra_dims ; if (use_extra_axis(m)) then - call write_field(IO_handle, more_axes(m)) + call IO_handle%write_field(more_axes(m)) endif ; enddo if (num_extra_dims > 0) then deallocate(use_extra_axis, more_axes) endif +end subroutine create_MOM_file + + +!> `reopen_MOM_file` wrapper for the legacy file handle, `file_type`. +!! NOTE: This function may be removed in a future release. +subroutine reopen_file(IO_handle, filename, vars, novars, fields, threading, & + timeunit, G, dG, GV, extra_axes, global_atts) + type(file_type), intent(inout) :: IO_handle + !< Handle for a file or fileset that is to be opened or reopened for + !! writing + character(len=*), intent(in) :: filename + !< full path to the file to create + type(vardesc), intent(in) :: vars(:) + !< structures describing fields written to filename + integer, intent(in) :: novars + !< number of fields written to filename + type(fieldtype), intent(inout) :: fields(:) + !< array of fieldtypes for each variable + integer, optional, intent(in) :: threading + !< SINGLE_FILE or MULTIPLE + real, optional, intent(in) :: timeunit + !< length of the units for time [s]. The default value is 86400.0, for 1 + !! day. + type(ocean_grid_type), optional, intent(in) :: G + !< ocean horizontal grid structure; G or dG is required if a 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 a new file + !! uses any horizontal grid axes. + type(verticalGrid_type), optional, intent(in) :: GV + !< ocean vertical grid structure, which is required if a new file uses any + !! vertical grid axes. + type(axis_info), optional, intent(in) :: extra_axes(:) + !< Types with information about some axes that might be used in this file + type(attribute_info), optional, intent(in) :: global_atts(:) + !< Global attributes to write to this file + + type(MOM_infra_file) :: mfile + !< Wrapper to MOM file + type(MOM_field), allocatable :: mfields(:) + !< Wrapper to MOM fields + integer :: i -end subroutine create_file + mfile%handle_infra = IO_handle + allocate(mfields(size(fields))) + + call reopen_MOM_file(mfile, filename, vars, novars, mfields, & + threading=threading, timeunit=timeunit, G=G, dG=dG, GV=GV, & + extra_axes=extra_axes, global_atts=global_atts) + + IO_handle = mfile%handle_infra + call get_file_fields(IO_handle, fields) +end subroutine reopen_file !> This routine opens an existing NetCDF file for output. If it !! does not find the file, a new file is created. It also sets up !! structures that describe this file and the variables that will !! later be written to this file. -subroutine reopen_file(IO_handle, filename, vars, novars, fields, threading, timeunit, & - G, dG, GV, extra_axes, global_atts) - type(file_type), intent(inout) :: IO_handle !< Handle for a file or fileset that is to be +subroutine reopen_MOM_file(IO_handle, filename, vars, novars, fields, & + threading, timeunit, G, dG, GV, extra_axes, global_atts) + class(MOM_file), intent(inout) :: IO_handle !< Handle for a file or fileset that is to be !! opened or reopened for writing character(len=*), intent(in) :: filename !< full path to the file to create type(vardesc), intent(in) :: vars(:) !< structures describing fields written to filename integer, intent(in) :: novars !< number of fields written to filename - type(fieldtype), intent(inout) :: fields(:) !< array of fieldtypes for each variable + type(MOM_field), intent(inout) :: fields(:) !< array of fieldtypes for each variable integer, optional, intent(in) :: threading !< SINGLE_FILE or MULTIPLE real, optional, intent(in) :: timeunit !< length of the units for time [s]. The !! default value is 86400.0, for 1 day. @@ -536,8 +657,9 @@ subroutine reopen_file(IO_handle, filename, vars, novars, fields, threading, tim inquire(file=check_name,EXIST=exists) if (.not.exists) then - call create_file(IO_handle, filename, vars, novars, fields, threading, timeunit, & - G=G, dG=dG, GV=GV, extra_axes=extra_axes, global_atts=global_atts) + call create_MOM_file(IO_handle, filename, vars, novars, fields, & + threading, timeunit, G=G, dG=dG, GV=GV, extra_axes=extra_axes, & + global_atts=global_atts) else domain_set = .false. @@ -551,40 +673,31 @@ subroutine reopen_file(IO_handle, filename, vars, novars, fields, threading, tim if (domain_set) one_file = (thread == SINGLE_FILE) if (one_file) then - call open_file(IO_handle, filename, APPEND_FILE, threading=thread) + call IO_handle%open(filename, APPEND_FILE, threading=thread) else - call open_file(IO_handle, filename, APPEND_FILE, MOM_domain=Domain) + call IO_handle%open(filename, APPEND_FILE, MOM_domain=Domain) endif - if (.not.file_is_open(IO_handle)) return + if (.not. IO_handle%file_is_open()) return - call get_file_info(IO_handle, nvar=nvar) + call IO_handle%get_file_info(nvar=nvar) if (nvar == -1) then write (mesg,*) "Reopening file ",trim(filename)," apparently had ",nvar,& " variables. Clobbering and creating file with ",novars," instead." call MOM_error(WARNING,"MOM_io: "//mesg) - call create_file(IO_handle, filename, vars, novars, fields, threading, timeunit, & - G=G, dG=dG, GV=GV, extra_axes=extra_axes, global_atts=global_atts) + call create_MOM_file(IO_handle, filename, vars, novars, fields, & + threading, timeunit, G=G, dG=dG, GV=GV, extra_axes=extra_axes, & + global_atts=global_atts) elseif (nvar /= novars) then write (mesg,*) "Reopening file ",trim(filename)," with ",novars,& " variables instead of ",nvar,"." call MOM_error(FATAL,"MOM_io: "//mesg) endif - if (nvar > 0) call get_file_fields(IO_handle, fields(1:nvar)) - - ! Check for inconsistent field names... -! do i=1,nvar -! call get_field_atts(fields(i), name) -! !if (trim(name) /= trim(vars%name)) then -! ! write (mesg, '("Reopening file ",a," variable ",a," is called ",a,".")',& -! ! trim(filename), trim(vars%name), trim(name)) -! ! call MOM_error(NOTE, "MOM_io: "//trim(mesg)) -! !endif -! enddo + if (nvar > 0) call IO_handle%get_file_fields(fields(1:nvar)) endif +end subroutine reopen_MOM_file -end subroutine reopen_file !> Return the index of sdtout if called from the root PE, or 0 for other PEs. integer function stdout_if_root() @@ -2089,9 +2202,8 @@ subroutine MOM_read_vector_3d(filename, u_fieldname, v_fieldname, u_data, v_data endif end subroutine MOM_read_vector_3d - !> 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, & +subroutine MOM_write_field_legacy_4d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, & fill_value, turns, scale) 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 @@ -2122,10 +2234,11 @@ subroutine MOM_write_field_4d(IO_handle, field_md, MOM_domain, field, tstamp, ti tile_count=tile_count, fill_value=fill_value) deallocate(field_rot) endif -end subroutine MOM_write_field_4d +end subroutine MOM_write_field_legacy_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, & +subroutine MOM_write_field_legacy_3d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, & fill_value, turns, scale) 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 @@ -2156,10 +2269,11 @@ subroutine MOM_write_field_3d(IO_handle, field_md, MOM_domain, field, tstamp, ti tile_count=tile_count, fill_value=fill_value) deallocate(field_rot) endif -end subroutine MOM_write_field_3d +end subroutine MOM_write_field_legacy_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, & +subroutine MOM_write_field_legacy_2d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, & fill_value, turns, scale) 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 @@ -2190,10 +2304,11 @@ subroutine MOM_write_field_2d(IO_handle, field_md, MOM_domain, field, tstamp, ti tile_count=tile_count, fill_value=fill_value) deallocate(field_rot) endif -end subroutine MOM_write_field_2d +end subroutine MOM_write_field_legacy_2d + !> Write a 1d field to an output file -subroutine MOM_write_field_1d(IO_handle, field_md, field, tstamp, fill_value, scale) +subroutine MOM_write_field_legacy_1d(IO_handle, field_md, field, tstamp, fill_value, scale) 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 @@ -2219,10 +2334,11 @@ subroutine MOM_write_field_1d(IO_handle, field_md, field, tstamp, fill_value, sc call write_field(IO_handle, field_md, array, tstamp=tstamp) deallocate(array) endif -end subroutine MOM_write_field_1d +end subroutine MOM_write_field_legacy_1d + !> Write a 0d field to an output file -subroutine MOM_write_field_0d(IO_handle, field_md, field, tstamp, fill_value, scale) +subroutine MOM_write_field_legacy_0d(IO_handle, field_md, field, tstamp, fill_value, scale) 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 @@ -2237,6 +2353,156 @@ subroutine MOM_write_field_0d(IO_handle, field_md, field, tstamp, fill_value, sc if (present(fill_value)) then ; if (field == fill_value) scaled_val = fill_value ; endif call write_field(IO_handle, field_md, scaled_val, tstamp=tstamp) +end subroutine MOM_write_field_legacy_0d + + +!> 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) + class(MOM_file), intent(inout) :: IO_handle !< Handle for a file that is open for writing + type(MOM_field), 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 + real, optional, intent(in) :: tstamp !< Model timestamp + integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) + real, optional, intent(in) :: fill_value !< Missing data fill value + integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data + real, optional, intent(in) :: scale !< A scaling factor that the field is + !! multiplied by before it is written + + real, allocatable :: field_rot(:,:,:,:) ! A rotated version of field, with the same units or rescaled + real :: scale_fac ! A scaling factor to use before writing the array + integer :: qturns ! The number of quarter turns through which to rotate field + + qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4) + scale_fac = 1.0 ; if (present(scale)) scale_fac = scale + + if ((qturns == 0) .and. (scale_fac == 1.0)) then + call IO_handle%write_field(field_md, MOM_domain, field, tstamp=tstamp, & + tile_count=tile_count, fill_value=fill_value) + else + call allocate_rotated_array(field, [1,1,1,1], qturns, field_rot) + call rotate_array(field, qturns, field_rot) + if (scale_fac /= 1.0) call rescale_comp_data(MOM_Domain, field_rot, scale_fac) + call IO_handle%write_field(field_md, MOM_domain, field_rot, tstamp=tstamp, & + tile_count=tile_count, fill_value=fill_value) + deallocate(field_rot) + endif +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) + class(MOM_file), intent(inout) :: IO_handle !< Handle for a file that is open for writing + type(MOM_field), 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 + real, optional, intent(in) :: tstamp !< Model timestamp + integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) + real, optional, intent(in) :: fill_value !< Missing data fill value + integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data + real, optional, intent(in) :: scale !< A scaling factor that the field is + !! multiplied by before it is written + + real, allocatable :: field_rot(:,:,:) ! A rotated version of field, with the same units or rescaled + real :: scale_fac ! A scaling factor to use before writing the array + integer :: qturns ! The number of quarter turns through which to rotate field + + qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4) + scale_fac = 1.0 ; if (present(scale)) scale_fac = scale + + if ((qturns == 0) .and. (scale_fac == 1.0)) then + call IO_handle%write_field(field_md, MOM_domain, field, tstamp=tstamp, & + tile_count=tile_count, fill_value=fill_value) + else + call allocate_rotated_array(field, [1,1,1], qturns, field_rot) + call rotate_array(field, qturns, field_rot) + if (scale_fac /= 1.0) call rescale_comp_data(MOM_Domain, field_rot, scale_fac) + call IO_handle%write_field(field_md, MOM_domain, field_rot, tstamp=tstamp, & + tile_count=tile_count, fill_value=fill_value) + deallocate(field_rot) + endif +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) + class(MOM_file), intent(inout) :: IO_handle !< Handle for a file that is open for writing + type(MOM_field), 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 + real, optional, intent(in) :: tstamp !< Model timestamp + integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1) + real, optional, intent(in) :: fill_value !< Missing data fill value + integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data + real, optional, intent(in) :: scale !< A scaling factor that the field is + !! multiplied by before it is written + + real, allocatable :: field_rot(:,:) ! A rotated version of field, with the same units + real :: scale_fac ! A scaling factor to use before writing the array + integer :: qturns ! The number of quarter turns through which to rotate field + + qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4) + scale_fac = 1.0 ; if (present(scale)) scale_fac = scale + + if ((qturns == 0) .and. (scale_fac == 1.0)) then + call IO_handle%write_field(field_md, MOM_domain, field, tstamp=tstamp, & + tile_count=tile_count, fill_value=fill_value) + else + call allocate_rotated_array(field, [1,1], qturns, field_rot) + call rotate_array(field, qturns, field_rot) + if (scale_fac /= 1.0) call rescale_comp_data(MOM_Domain, field_rot, scale_fac) + call IO_handle%write_field(field_md, MOM_domain, field_rot, tstamp=tstamp, & + tile_count=tile_count, fill_value=fill_value) + deallocate(field_rot) + endif +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) + class(MOM_file), intent(inout) :: IO_handle !< Handle for a file that is open for writing + type(MOM_field), 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) :: fill_value !< Missing data fill value + real, optional, intent(in) :: scale !< A scaling factor that the field is + !! multiplied by before it is written + + real, dimension(:), allocatable :: array ! A rescaled copy of field + real :: scale_fac ! A scaling factor to use before writing the array + integer :: i + + scale_fac = 1.0 ; if (present(scale)) scale_fac = scale + + if (scale_fac == 1.0) then + call IO_handle%write_field(field_md, field, tstamp=tstamp) + else + allocate(array(size(field))) + array(:) = scale_fac * field(:) + if (present(fill_value)) then + do i=1,size(field) ; if (field(i) == fill_value) array(i) = fill_value ; enddo + endif + call IO_handle%write_field(field_md, array, tstamp=tstamp) + deallocate(array) + endif +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) + class(MOM_file), intent(inout) :: IO_handle !< Handle for a file that is open for writing + type(MOM_field), 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) :: fill_value !< Missing data fill value + real, optional, intent(in) :: scale !< A scaling factor that the field is + !! multiplied by before it is written + real :: scaled_val ! A rescaled copy of field + + scaled_val = field + if (present(scale)) scaled_val = scale*field + if (present(fill_value)) then ; if (field == fill_value) scaled_val = fill_value ; endif + + call IO_handle%write_field(field_md, scaled_val, tstamp=tstamp) end subroutine MOM_write_field_0d !> Given filename and fieldname, this subroutine returns the size of the field in the file diff --git a/src/framework/MOM_io_file.F90 b/src/framework/MOM_io_file.F90 new file mode 100644 index 0000000000..68c0f33f07 --- /dev/null +++ b/src/framework/MOM_io_file.F90 @@ -0,0 +1,1654 @@ +!> This module contains the MOM file handler types +module MOM_io_file + +! This file is part of MOM6. See LICENSE.md for the license. + +use, intrinsic :: iso_fortran_env, only : int64 + +use MOM_domains, only : MOM_domain_type, domain1D +use MOM_io_infra, only : file_type, get_file_info, get_file_fields +use MOM_io_infra, only : open_file, close_file, flush_file +use MOM_io_infra, only : fms2_file_is_open => file_is_open +use MOM_io_infra, only : fieldtype +use MOM_io_infra, only : get_file_times, axistype +use MOM_io_infra, only : write_field, write_metadata +use MOM_io_infra, only : get_field_atts +use MOM_io_infra, only : read_field_chksum + +use MOM_netcdf, only : netcdf_file_type +use MOM_netcdf, only : netcdf_axis +use MOM_netcdf, only : netcdf_field +use MOM_netcdf, only : open_netcdf_file +use MOM_netcdf, only : close_netcdf_file +use MOM_netcdf, only : flush_netcdf_file +use MOM_netcdf, only : register_netcdf_axis +use MOM_netcdf, only : register_netcdf_field +use MOM_netcdf, only : write_netcdf_field +use MOM_netcdf, only : write_netcdf_axis +use MOM_netcdf, only : write_netcdf_attribute +use MOM_netcdf, only : get_netcdf_size +use MOM_netcdf, only : get_netcdf_fields + +use MOM_error_handler, only : MOM_error, FATAL +use MOM_error_handler, only : is_root_PE + +implicit none ; private + +public :: MOM_file +public :: MOM_infra_file +public :: MOM_netcdf_file +public :: MOM_axis +public :: MOM_field + + +! Internal types + +! NOTE: MOM_axis and MOM_field do not represent the actual axes and +! fields stored in the file. They are only very thin wrappers to the keys (as +! strings) used to reference the associated object inside the MOM_file. + +!> Handle for axis in MOM file +type :: MOM_axis + character(len=:), allocatable :: label + !< Identifier for the axis in handle's list +end type MOM_axis + + +!> Linked list of framework axes +type :: axis_list_infra + private + type(axis_node_infra), pointer :: head => null() + !< Head of axis linked list + type(axis_node_infra), pointer :: tail => null() + !< Tail of axis linked list +contains + !> Initialize the framework axis list + procedure :: init => initialize_axis_list_infra + !> Append a new axis to the framework axis list + procedure :: append => append_axis_list_infra + !> Get an axis from the framework axis list + procedure :: get => get_axis_list_infra + !> Deallocate the framework axis list + procedure :: finalize => finalize_axis_list_infra +end type axis_list_infra + + +!> Framework axis linked list node +type :: axis_node_infra + private + character(len=:), allocatable :: label + !< Axis identifier + type(axis_node_infra), pointer :: next => null() + !< Pointer to next axis node + type(axistype) :: axis + !< Axis node contents +end type axis_node_infra + + +!> Linked list of framework axes +type :: axis_list_nc + private + type(axis_node_nc), pointer :: head => null() + !< Head of axis linked list + type(axis_node_nc), pointer :: tail => null() + !< Tail of axis linked list +contains + !> Initialize the netCDF axis list + procedure :: init => initialize_axis_list_nc + !> Append a new axis to the netCDF axis list + procedure :: append => append_axis_list_nc + !> Get an axis from the netCDF axis list + procedure :: get => get_axis_list_nc + !> Deallocate the netCDF axis list + procedure :: finalize => finalize_axis_list_nc +end type axis_list_nc + + +!> Framework axis linked list node +type :: axis_node_nc + private + character(len=:), allocatable :: label + !< Axis identifier + type(axis_node_nc), pointer :: next => null() + !< Pointer to next axis node + type(netcdf_axis) :: axis + !< Axis node contents +end type axis_node_nc + + +!> Handle for field in MOM file +type :: MOM_field + character(len=:), allocatable :: label + !< Identifier for the field in the handle's list +end type MOM_field + + +!> Linked list of framework fields +type :: field_list_infra + private + type(field_node_infra), pointer :: head => null() + !< Head of field linked list + type(field_node_infra), pointer :: tail => null() + !< Tail of field linked list +contains + !> Initialize the framework field list + procedure :: init => initialize_field_list_infra + !> Append a new axis to the framework field list + procedure :: append => append_field_list_infra + !> Get an axis from the framework field list + procedure :: get => get_field_list_infra + !> Deallocate the framework field list + procedure :: finalize => finalize_field_list_infra +end type field_list_infra + + +!> Framework field linked list node +type :: field_node_infra + private + character(len=:), allocatable :: label + !< Field identifier + type(fieldtype) :: field + !< Field node contents + type(field_node_infra), pointer :: next => null() + !< Pointer to next field node +end type field_node_infra + + +!> Linked list of framework fields +type :: field_list_nc + private + type(field_node_nc), pointer :: head => null() + !< Head of field linked list + type(field_node_nc), pointer :: tail => null() + !< Tail of field linked list +contains + !> Initialize the netCDF field list + procedure :: init => initialize_field_list_nc + !> Append a new axis to the netCDF field list + procedure :: append => append_field_list_nc + !> Get an axis from the netCDF field list + procedure :: get => get_field_list_nc + !> Deallocate the netCDF field list + procedure :: finalize => finalize_field_list_nc +end type field_list_nc + + +!> Framework field linked list node +type :: field_node_nc + private + character(len=:), allocatable :: label + !< Field identifier + type(netcdf_field) :: field + !< Field node contents + type(field_node_nc), pointer :: next => null() + !< Pointer to next field node +end type field_node_nc + + +!> Generic MOM file abstraction for common operations +type, abstract :: MOM_file + private + + contains + + !> Open a file and connect to the MOM_file object + procedure(i_open_file), deferred :: open + !> Close the MOM file + procedure(i_close_file), deferred :: close + !> Flush buffered output to the MOM file + procedure(i_flush_file), deferred :: flush + + !> Register an axis to the MOM file + procedure(i_register_axis), deferred :: register_axis + !> Register a field to the MOM file + procedure(i_register_field), deferred :: register_field + !> Write metadata to the MOM file + procedure(i_write_attribute), deferred :: write_attribute + + !> Write field to a MOM file + generic :: write_field => & + write_field_4d, & + write_field_3d, & + write_field_2d, & + write_field_1d, & + write_field_0d, & + write_field_axis + + !> Write a 4D field to the MOM file + procedure(i_write_field_4d), deferred :: write_field_4d + !> Write a 3D field to the MOM file + procedure(i_write_field_3d), deferred :: write_field_3d + !> Write a 2D field to the MOM file + procedure(i_write_field_2d), deferred :: write_field_2d + !> Write a 1D field to the MOM file + procedure(i_write_field_1d), deferred :: write_field_1d + !> Write a 0D field to the MOM file + procedure(i_write_field_0d), deferred :: write_field_0d + !> Write an axis field to the MOM file + procedure(i_write_field_axis), deferred :: write_field_axis + + !> Return true if MOM file has been opened + procedure(i_file_is_open), deferred :: file_is_open + !> Return number of dimensions, variables, or time levels in a MOM file + procedure(i_get_file_info), deferred :: get_file_info + !> Get field objects from a MOM file + procedure(i_get_file_fields), deferred :: get_file_fields + !> Get attributes from a field + procedure(i_get_field_atts), deferred :: get_field_atts + !> Get checksum from a field + procedure(i_read_field_chksum), deferred :: read_field_chksum +end type MOM_file + + +!> MOM file from the supporting framework ("infra") layer +type, extends(MOM_file) :: MOM_infra_file + private + + ! NOTE: This will be made private after the API transition + type(file_type), public :: handle_infra + !< Framework-specific file handler content + type(axis_list_infra) :: axes + !< List of axes in file + type(field_list_infra) :: fields + !< List of fields in file + + contains + + !> Open a framework file and connect to the MOM_file object + procedure :: open => open_file_infra + !> Close the MOM framework file + procedure :: close => close_file_infra + !> Flush buffered output to the MOM framework file + procedure :: flush => flush_file_infra + + !> Register an axis to the MOM framework file + procedure :: register_axis => register_axis_infra + !> Register a field to the MOM framework file + procedure :: register_field => register_field_infra + !> Write global metadata to the MOM framework file + procedure :: write_attribute => write_attribute_infra + + !> Write a 4D field to the MOM framework file + procedure :: write_field_4d => write_field_4d_infra + !> Write a 3D field to the MOM framework file + procedure :: write_field_3d => write_field_3d_infra + !> Write a 2D field to the MOM framework file + procedure :: write_field_2d => write_field_2d_infra + !> Write a 1D field to the MOM framework file + procedure :: write_field_1d => write_field_1d_infra + !> Write a 0D field to the MOM framework file + procedure :: write_field_0d => write_field_0d_infra + !> Write an axis field to the MOM framework file + procedure :: write_field_axis => write_field_axis_infra + + !> Return true if MOM infra file has been opened + procedure :: file_is_open => file_is_open_infra + !> Return number of dimensions, variables, or time levels in a MOM infra file + procedure :: get_file_info => get_file_info_infra + !> Get field metadata from a MOM infra file + procedure :: get_file_fields => get_file_fields_infra + !> Get attributes from a field + procedure :: get_field_atts => get_field_atts_infra + !> Get checksum from a field + procedure :: read_field_chksum => read_field_chksum_infra + + ! MOM_infra_file methods + ! NOTE: These could naturally reside in MOM_file but is currently not needed. + + !> Get time levels of a MOM framework file + procedure :: get_file_times => get_file_times_infra + + !> Get the fields as fieldtypes from a file + procedure :: get_file_fieldtypes + ! NOTE: This is provided to support the legacy API and may be removed. +end type MOM_infra_file + + +!> MOM file using netCDF backend +type, extends(MOM_file) :: MOM_netcdf_file + private + + !> Framework-specific file handler content + type(netcdf_file_type) :: handle_nc + !> List of netCDF axes + type(axis_list_nc) :: axes + !> List of netCDF fields + type(field_list_nc) :: fields + !> True if the file has been opened + logical :: is_open = .false. + + contains + + !> Open a framework file and connect to the MOM_netcdf_file object + procedure :: open => open_file_nc + !> Close the MOM netcdf file + procedure :: close => close_file_nc + !> Flush buffered output to the MOM netcdf file + procedure :: flush => flush_file_nc + + !> Register an axis to the MOM netcdf file + procedure :: register_axis => register_axis_nc + !> Register a field to the MOM netcdf file + procedure :: register_field => register_field_nc + !> Write global metadata to the MOM netcdf file + procedure :: write_attribute => write_attribute_nc + + !> Write a 4D field to the MOM netcdf file + procedure :: write_field_4d => write_field_4d_nc + !> Write a 3D field to the MOM netcdf file + procedure :: write_field_3d => write_field_3d_nc + !> Write a 2D field to the MOM netcdf file + procedure :: write_field_2d => write_field_2d_nc + !> Write a 1D field to the MOM netcdf file + procedure :: write_field_1d => write_field_1d_nc + !> Write a 0D field to the MOM netcdf file + procedure :: write_field_0d => write_field_0d_nc + !> Write an axis field to the MOM netcdf file + procedure :: write_field_axis => write_field_axis_nc + + !> Return true if MOM netcdf file has been opened + procedure :: file_is_open => file_is_open_nc + !> Return number of dimensions, variables, or time levels in a MOM netcdf file + procedure :: get_file_info => get_file_info_nc + !> Get field metadata from a MOM netcdf file + procedure :: get_file_fields => get_file_fields_nc + !> Get attributes from a netCDF field + procedure :: get_field_atts => get_field_atts_nc + !> Get checksum from a netCDF field + procedure :: read_field_chksum => read_field_chksum_nc +end type MOM_netcdf_file + + +interface + !> Interface for opening a MOM file + subroutine i_open_file(handle, filename, action, MOM_domain, threading, fileset) + import :: MOM_file, MOM_domain_type + + class(MOM_file), intent(inout) :: handle + !< The handle for the opened file + character(len=*), intent(in) :: filename + !< The path name of the file being opened + integer, optional, intent(in) :: action + !< A flag indicating whether the file can be read or written to and how + !! to handle existing files. The default is WRITE_ONLY. + type(MOM_domain_type), optional, intent(in) :: MOM_Domain + !< A MOM_Domain that describes the decomposition + integer, optional, intent(in) :: threading + !< A flag indicating whether one (SINGLE_FILE) or multiple PEs (MULTIPLE) + !! participate in I/O. With the default, the root PE does I/O. + integer, optional, intent(in) :: fileset + !< A flag indicating whether multiple PEs doing I/O due to + !! threading=MULTIPLE write to the same file (SINGLE_FILE) or to one file + !! per PE (MULTIPLE, the default). + end subroutine i_open_file + + + !> Interface for closing a MOM file + subroutine i_close_file(handle) + import :: MOM_file + class(MOM_file), intent(inout) :: handle + !< The MOM file to be closed + end subroutine i_close_file + + + !> Interface for flushing I/O in a MOM file + subroutine i_flush_file(handle) + import :: MOM_file + class(MOM_file), intent(in) :: handle + !< The MOM file to be flushed + end subroutine i_flush_file + + + !> Interface to register an axis to a MOM file + function i_register_axis(handle, label, units, longname, cartesian, sense, & + domain, data, edge_axis, calendar) result(axis) + import :: MOM_file, MOM_axis, domain1D + + class(MOM_file), intent(inout) :: handle + !< Handle for a file that is open for writing + character(len=*), intent(in) :: label + !< The name in the file of this axis + character(len=*), intent(in) :: units + !< The units of this axis + character(len=*), intent(in) :: longname + !< The long description of this axis + character(len=*), optional, intent(in) :: cartesian + !< A variable indicating which direction this axis corresponds with. + !! Valid values include 'X', 'Y', 'Z', 'T', and 'N' for none. + integer, optional, intent(in) :: sense + !< This is 1 for axes whose values increase upward, or -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 + type(MOM_axis) :: axis + !< IO handle for axis in MOM_file + end function i_register_axis + + + !> Interface to register a field to a netCDF file + function i_register_field(handle, axes, label, units, longname, & + pack, standard_name, checksum) result(field) + import :: MOM_file, MOM_axis, MOM_field, int64 + class(MOM_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_axis), intent(in) :: axes(:) + !< Handles for the axis used for this variable + character(len=*), intent(in) :: label + !< 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 + integer, optional, intent(in) :: pack + !< A precision reduction factor with which the variable. The default, 1, + !! has no reduction, but 2 is not uncommon. + character(len=*), optional, intent(in) :: standard_name + !< The standard (e.g., CMOR) name for this variable + integer(kind=int64), dimension(:), optional, intent(in) :: checksum + !< Checksum values that can be used to verify reads. + type(MOM_field) :: field + !< IO handle for field in MOM_file + end function i_register_field + + + !> Interface for writing global metata to a MOM file + subroutine i_write_attribute(handle, name, attribute) + import :: MOM_file + class(MOM_file), intent(in) :: handle + !< Handle for a file that is open for writing + character(len=*), intent(in) :: name + !< The name in the file of this global attribute + character(len=*), intent(in) :: attribute + !< The value of this attribute + end subroutine i_write_attribute + + + !> Interface to write_field_4d() + subroutine i_write_field_4d(handle, field_md, MOM_domain, field, tstamp, & + tile_count, fill_value) + import :: MOM_file, MOM_field, MOM_domain_type + class(MOM_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_field), intent(in) :: field_md + !< Field type with metadata + type(MOM_domain_type), intent(in) :: MOM_domain + !< The MOM_Domain that describes the decomposition + real, intent(inout) :: field(:,:,:,:) + !< Field to write + 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 + end subroutine i_write_field_4d + + + !> Interface to write_field_3d() + subroutine i_write_field_3d(handle, field_md, MOM_domain, field, tstamp, & + tile_count, fill_value) + import :: MOM_file, MOM_field, MOM_domain_type + class(MOM_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_field), intent(in) :: field_md + !< Field type with metadata + type(MOM_domain_type), intent(in) :: MOM_domain + !< The MOM_Domain that describes the decomposition + real, intent(inout) :: field(:,:,:) + !< Field to write + 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 + end subroutine i_write_field_3d + + + !> Interface to write_field_2d() + subroutine i_write_field_2d(handle, field_md, MOM_domain, field, tstamp, & + tile_count, fill_value) + import :: MOM_file, MOM_field, MOM_domain_type + class(MOM_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_field), 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 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 + end subroutine i_write_field_2d + + + !> Interface to write_field_1d() + subroutine i_write_field_1d(handle, field_md, field, tstamp) + import :: MOM_file, MOM_field + class(MOM_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_field), 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 + end subroutine i_write_field_1d + + + !> Interface to write_field_0d() + subroutine i_write_field_0d(handle, field_md, field, tstamp) + import :: MOM_file, MOM_field + class(MOM_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_field), 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 + end subroutine i_write_field_0d + + + !> Interface to write_field_axis() + subroutine i_write_field_axis(handle, axis) + import :: MOM_file, MOM_axis + class(MOM_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_axis), intent(in) :: axis + !< An axis type variable with information to write + end subroutine i_write_field_axis + + + !> Interface to file_is_open() + logical function i_file_is_open(handle) + import :: MOM_file + class(MOM_file), intent(in) :: handle + !< Handle to a file to inquire about + end function i_file_is_open + + + !> Interface to get_file_info() + subroutine i_get_file_info(handle, ndim, nvar, ntime) + import :: MOM_file + class(MOM_file), intent(in) :: 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) :: ntime + !< The number of time levels in the file + end subroutine i_get_file_info + + + !> Interface to get_file_fields() + subroutine i_get_file_fields(handle, fields) + import :: MOM_file, MOM_field + class(MOM_file), intent(inout) :: handle + !< Handle for a file that is open for I/O + type(MOM_field), dimension(:), intent(inout) :: fields + !< Field-type descriptions of all of the variables in a file. + end subroutine i_get_file_fields + + + !> Interface to get_field_atts() + subroutine i_get_field_atts(handle, field, name, units, longname, checksum) + import :: MOM_file, MOM_field, int64 + class(MOM_file), intent(in) :: handle + !< File where field is stored + type(MOM_field), intent(in) :: field + !< The field to extract information from + character(len=*), optional, intent(out) :: name + !< The variable name + character(len=*), optional, intent(out) :: units + !< The units of the variable + character(len=*), optional, intent(out) :: longname + !< The long name of the variable + integer(kind=int64), optional, intent(out) :: checksum(:) + !< The checksums of the variable in a file + end subroutine i_get_field_atts + + + !> Interface to read_field_chksum + subroutine i_read_field_chksum(handle, field, chksum, valid_chksum) + import :: MOM_file, MOM_field, int64 + class(MOM_file), intent(in) :: handle + !< File where field is stored + type(MOM_field), 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 + end subroutine i_read_field_chksum +end interface + +contains + +!> Initialize the linked list of framework axes +subroutine initialize_axis_list_infra(list) + class(axis_list_infra), intent(inout) :: list + + ! Pre-allocate the first node and set the tail to this empty node + allocate(list%head) + list%tail => list%head +end subroutine initialize_axis_list_infra + + +!> Append a new axis to the list +subroutine append_axis_list_infra(list, axis, label) + class(axis_list_infra), intent(inout) :: list + type(axistype), intent(in) :: axis + character(len=*), intent(in) :: label + + type(axis_node_infra), pointer :: empty_node + + ! Transfer value to tail + list%tail%label = label + list%tail%axis = axis + + ! Extend list to next empty node + allocate(empty_node) + list%tail%next => empty_node + list%tail => empty_node +end subroutine append_axis_list_infra + + +!> Get axis based on label +function get_axis_list_infra(list, label) result(axis) + class(axis_list_infra), intent(in) :: list + character(len=*), intent(in) :: label + type(axistype) :: axis + + type(axis_node_infra), pointer :: node + + ! NOTE: The tail is a pre-allocated empty node, so we check node%next + node => list%head + do while(associated(node%next)) + if (node%label == label) exit + node => node%next + enddo + if (.not. associated(node)) & + call MOM_error(FATAL, "axis associated with " // label // " not found.") + + axis = node%axis +end function get_axis_list_infra + + +!> Deallocate axes of list +subroutine finalize_axis_list_infra(list) + class(axis_list_infra), intent(inout) :: list + + type(axis_node_infra), pointer :: node, next_node + + node => list%head + do while(associated(node)) + next_node => node + node => node%next + deallocate(next_node) + enddo +end subroutine finalize_axis_list_infra + + +!> Initialize the linked list of framework axes +subroutine initialize_axis_list_nc(list) + class(axis_list_nc), intent(inout) :: list + + ! Pre-allocate the first node and set the tail to this empty node + allocate(list%head) + list%tail => list%head +end subroutine initialize_axis_list_nc + + +!> Append a new axis to the list +subroutine append_axis_list_nc(list, axis, label) + class(axis_list_nc), intent(inout) :: list + type(netcdf_axis), intent(in) :: axis + character(len=*), intent(in) :: label + + type(axis_node_nc), pointer :: empty_node + + ! Transfer value to tail + list%tail%label = label + list%tail%axis = axis + + ! Extend list to next empty node + allocate(empty_node) + list%tail%next => empty_node + list%tail => empty_node +end subroutine append_axis_list_nc + + +!> Get axis based on label +function get_axis_list_nc(list, label) result(axis) + class(axis_list_nc), intent(in) :: list + character(len=*), intent(in) :: label + type(netcdf_axis) :: axis + + type(axis_node_nc), pointer :: node + + ! NOTE: The tail is a pre-allocated empty node, so we check node%next + node => list%head + do while(associated(node%next)) + if (node%label == label) exit + node => node%next + enddo + if (.not. associated(node)) & + call MOM_error(FATAL, "axis associated with " // label // " not found.") + + axis = node%axis +end function get_axis_list_nc + + +!> Deallocate axes of list +subroutine finalize_axis_list_nc(list) + class(axis_list_nc), intent(inout) :: list + + type(axis_node_nc), pointer :: node, next_node + + node => list%head + do while(associated(node)) + next_node => node + node => node%next + deallocate(next_node) + enddo +end subroutine finalize_axis_list_nc + + +!> Initialize the linked list of framework axes +subroutine initialize_field_list_infra(list) + class(field_list_infra), intent(inout) :: list + + ! Pre-allocate the first node and set the tail to this empty node + allocate(list%head) + list%tail => list%head +end subroutine initialize_field_list_infra + + +!> Append a new field to the list +subroutine append_field_list_infra(list, field, label) + class(field_list_infra), intent(inout) :: list + type(fieldtype), intent(in) :: field + character(len=*), intent(in) :: label + + type(field_node_infra), pointer :: empty_node + + ! Transfer value to tail + list%tail%label = label + list%tail%field = field + + ! Extend list to next empty node + allocate(empty_node) + list%tail%next => empty_node + list%tail => empty_node +end subroutine append_field_list_infra + + +!> Get axis based on label +function get_field_list_infra(list, label) result(field) + class(field_list_infra), intent(in) :: list + character(len=*), intent(in) :: label + type(fieldtype) :: field + + type(field_node_infra), pointer :: node + + ! NOTE: The tail is a pre-allocated empty node, so we check node%next + node => list%head + do while(associated(node%next)) + if (node%label == label) exit + node => node%next + enddo + if (.not. associated(node)) & + call MOM_error(FATAL, "field associated with " // label // " not found.") + + field = node%field +end function get_field_list_infra + + +!> Deallocate fields of list +subroutine finalize_field_list_infra(list) + class(field_list_infra), intent(inout) :: list + + type(field_node_infra), pointer :: node, next_node + + node => list%head + do while(associated(node)) + next_node => node + node => node%next + deallocate(next_node) + enddo +end subroutine finalize_field_list_infra + + +!> Initialize the linked list of framework axes +subroutine initialize_field_list_nc(list) + class(field_list_nc), intent(inout) :: list + + ! Pre-allocate the first node and set the tail to this empty node + allocate(list%head) + list%tail => list%head +end subroutine initialize_field_list_nc + + +!> Append a new field to the list +subroutine append_field_list_nc(list, field, label) + class(field_list_nc), intent(inout) :: list + type(netcdf_field), intent(in) :: field + character(len=*), intent(in) :: label + + type(field_node_nc), pointer :: empty_node + + ! Transfer value to tail + list%tail%label = label + list%tail%field = field + + ! Extend list to next empty node + allocate(empty_node) + list%tail%next => empty_node + list%tail => empty_node +end subroutine append_field_list_nc + + +!> Get axis based on label +function get_field_list_nc(list, label) result(field) + class(field_list_nc), intent(in) :: list + character(len=*), intent(in) :: label + type(netcdf_field) :: field + + type(field_node_nc), pointer :: node + + ! NOTE: The tail is a pre-allocated empty node, so we check node%next + node => list%head + do while(associated(node%next)) + if (node%label == label) exit + node => node%next + enddo + if (.not. associated(node)) & + call MOM_error(FATAL, "field associated with " // label // " not found.") + + field = node%field +end function get_field_list_nc + + +!> Deallocate fields of list +subroutine finalize_field_list_nc(list) + class(field_list_nc), intent(inout) :: list + + type(field_node_nc), pointer :: node, next_node + + node => list%head + do while(associated(node)) + next_node => node + node => node%next + deallocate(next_node) + enddo +end subroutine finalize_field_list_nc + + +!> Open a MOM framework file +subroutine open_file_infra(handle, filename, action, MOM_domain, threading, fileset) + class(MOM_infra_file), intent(inout) :: handle + character(len=*), intent(in) :: filename + integer, intent(in), optional :: action + type(MOM_domain_type), optional, intent(in) :: MOM_domain + integer, intent(in), optional :: threading + integer, intent(in), optional :: fileset + + call open_file(handle%handle_infra, filename, action=action, & + MOM_domain=MOM_domain, threading=threading, fileset=fileset) + + call handle%axes%init() + call handle%fields%init() +end subroutine open_file_infra + +!> Close a MOM framework file +subroutine close_file_infra(handle) + class(MOM_infra_file), intent(inout) :: handle + + call close_file(handle%handle_infra) + call handle%axes%finalize() + call handle%fields%finalize() +end subroutine close_file_infra + +!> Flush the buffer of a MOM framework file +subroutine flush_file_infra(handle) + class(MOM_infra_file), intent(in) :: handle + + call flush_file(handle%handle_infra) +end subroutine flush_file_infra + + +!> Register an axis to the MOM framework file +function register_axis_infra(handle, label, units, longname, & + cartesian, sense, domain, data, edge_axis, calendar) result(axis) + + class(MOM_infra_file), intent(inout) :: handle + !< Handle for a file that is open for writing + character(len=*), intent(in) :: label + !< The name in the file of this axis + character(len=*), intent(in) :: units + !< The units of this axis + character(len=*), intent(in) :: longname + !< The long description of this axis + character(len=*), optional, intent(in) :: cartesian + !< A variable indicating which direction this axis corresponds with. + !! Valid values include 'X', 'Y', 'Z', 'T', and 'N' for none. + integer, optional, intent(in) :: sense + !< This is 1 for axes whose values increase upward, or -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 + type(MOM_axis) :: axis + !< The axis type where this information is stored + + type(axistype) :: ax_infra + + ! Create new infra axis and assign to pre-allocated tail of axes + call write_metadata(handle%handle_infra, ax_infra, label, units, longname, & + cartesian=cartesian, sense=sense, domain=domain, data=data, & + edge_axis=edge_axis, calendar=calendar) + + call handle%axes%append(ax_infra, label) + axis%label = label +end function register_axis_infra + + +!> Register a field to the MOM framework file +function register_field_infra(handle, axes, label, units, longname, pack, & + standard_name, checksum) result(field) + class(MOM_infra_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_axis), dimension(:), intent(in) :: axes + !< Handles for the axis used for this variable + character(len=*), intent(in) :: label + !< 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 + integer, optional, intent(in) :: pack + !< A precision reduction factor with which the variable. The default, 1, + !! has no reduction, but 2 is not uncommon. + character(len=*), optional, intent(in) :: standard_name + !< The standard (e.g., CMOR) name for this variable + integer(kind=int64), dimension(:), optional, intent(in) :: checksum + !< Checksum values that can be used to verify reads. + type(MOM_field) :: field + !< The field type where this information is stored + + type(fieldtype) :: field_infra + type(axistype), allocatable :: field_axes(:) + integer :: i + + ! Construct array of framework axes + allocate(field_axes(size(axes))) + do i = 1, size(axes) + field_axes(i) = handle%axes%get(axes(i)%label) + enddo + + call write_metadata(handle%handle_infra, field_infra, field_axes, label, & + units, longname, pack=pack, standard_name=standard_name, checksum=checksum) + + call handle%fields%append(field_infra, label) + field%label = label +end function register_field_infra + + +!> Write a 4D field to the MOM framework file +subroutine write_field_4d_infra(handle, field_md, MOM_domain, field, tstamp, & + tile_count, fill_value) + class(MOM_infra_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_field), intent(in) :: field_md + !< Field type with metadata + type(MOM_domain_type), intent(in) :: MOM_domain + !< The MOM_Domain that describes the decomposition + real, intent(inout) :: field(:,:,:,:) + !< Field to write + 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 + + type(fieldtype) :: field_infra + + field_infra = handle%fields%get(field_md%label) + call write_field(handle%handle_infra, field_infra, MOM_domain, field, & + tstamp=tstamp, tile_count=tile_count, fill_value=fill_value) +end subroutine write_field_4d_infra + + +!> Write a 3D field to the MOM framework file +subroutine write_field_3d_infra(handle, field_md, MOM_domain, field, tstamp, & + tile_count, fill_value) + class(MOM_infra_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_field), intent(in) :: field_md + !< Field type with metadata + type(MOM_domain_type), intent(in) :: MOM_domain + !< The MOM_Domain that describes the decomposition + real, intent(inout) :: field(:,:,:) + !< Field to write + 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 + + type(fieldtype) :: field_infra + + field_infra = handle%fields%get(field_md%label) + call write_field(handle%handle_infra, field_infra, MOM_domain, field, & + tstamp=tstamp, tile_count=tile_count, fill_value=fill_value) +end subroutine write_field_3d_infra + + +!> Write a 2D field to the MOM framework file +subroutine write_field_2d_infra(handle, field_md, MOM_domain, field, tstamp, & + tile_count, fill_value) + class(MOM_infra_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_field), 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 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 + + type(fieldtype) :: field_infra + + field_infra = handle%fields%get(field_md%label) + call write_field(handle%handle_infra, field_infra, MOM_domain, field, & + tstamp=tstamp, tile_count=tile_count, fill_value=fill_value) +end subroutine write_field_2d_infra + + +!> Write a 1D field to the MOM framework file +subroutine write_field_1d_infra(handle, field_md, field, tstamp) + class(MOM_infra_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_field), 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 + + type(fieldtype) :: field_infra + + field_infra = handle%fields%get(field_md%label) + call write_field(handle%handle_infra, field_infra, field, tstamp=tstamp) +end subroutine write_field_1d_infra + + +!> Write a 0D field to the MOM framework file +subroutine write_field_0d_infra(handle, field_md, field, tstamp) + class(MOM_infra_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_field), 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 + + type(fieldtype) :: field_infra + + field_infra = handle%fields%get(field_md%label) + call write_field(handle%handle_infra, field_infra, field, tstamp=tstamp) +end subroutine write_field_0d_infra + + +!> Write an axis field to the MOM framework file +subroutine write_field_axis_infra(handle, axis) + class(MOM_infra_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_axis), intent(in) :: axis + !< An axis type variable with information to write + + type(axistype) :: axis_infra + !< An axis type variable with information to write + + axis_infra = handle%axes%get(axis%label) + call write_field(handle%handle_infra, axis_infra) +end subroutine write_field_axis_infra + + +!> Write global metadata to the MOM framework file +subroutine write_attribute_infra(handle, name, attribute) + class(MOM_infra_file), intent(in) :: handle + !< Handle for a file that is open for writing + character(len=*), intent(in) :: name + !< The name in the file of this global attribute + character(len=*), intent(in) :: attribute + !< The value of this attribute + + call write_metadata(handle%handle_infra, name, attribute) +end subroutine write_attribute_infra + + +!> True if the framework file has been opened +logical function file_is_open_infra(handle) + class(MOM_infra_file), intent(in) :: handle + !< Handle to a file to inquire about + + file_is_open_infra = fms2_file_is_open(handle%handle_infra) +end function file_is_open_infra + + +!> Return number of dimensions, variables, or time levels in a MOM infra file +subroutine get_file_info_infra(handle, ndim, nvar, ntime) + class(MOM_infra_file), intent(in) :: 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) :: ntime + !< The number of time levels in the file + + call get_file_info(handle%handle_infra, ndim, nvar, ntime) +end subroutine get_file_info_infra + + +!> Return the field metadata associated with a MOM framework file +subroutine get_file_fields_infra(handle, fields) + class(MOM_infra_file), intent(inout) :: handle + !< Handle for a file that is open for I/O + type(MOM_field), intent(inout) :: fields(:) + !< Field-type descriptions of all of the variables in a file. + + type(fieldtype), allocatable :: fields_infra(:) + integer :: i + character(len=64) :: label + + allocate(fields_infra(size(fields))) + call get_file_fields(handle%handle_infra, fields_infra) + + do i = 1, size(fields) + call get_field_atts(fields_infra(i), name=label) + call handle%fields%append(fields_infra(i), trim(label)) + fields(i)%label = trim(label) + enddo +end subroutine get_file_fields_infra + + +!> Get time levels of a MOM framework file +subroutine get_file_times_infra(handle, time_values, ntime) + class(MOM_infra_file), intent(in) :: 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 + + call get_file_times(handle%handle_infra, time_values, ntime=ntime) +end subroutine get_file_times_infra + + +!> Get attributes from a field +subroutine get_field_atts_infra(handle, field, name, units, longname, checksum) + class(MOM_infra_file), intent(in) :: handle + !< File where field is stored + type(MOM_field), intent(in) :: field + !< The field to extract information from + character(len=*), optional, intent(out) :: name + !< The variable name + character(len=*), optional, intent(out) :: units + !< The units of the variable + character(len=*), optional, intent(out) :: longname + !< The long name of the variable + integer(kind=int64), optional, intent(out) :: checksum(:) + !< The checksums of the variable in a file + + type(fieldtype) :: field_infra + + field_infra = handle%fields%get(field%label) + call get_field_atts(field_infra, name, units, longname, checksum) +end subroutine get_field_atts_infra + + +!> Interface to read_field_chksum +subroutine read_field_chksum_infra(handle, field, chksum, valid_chksum) + class(MOM_infra_file), intent(in) :: handle + !< File where field is stored + type(MOM_field), 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 + + type(fieldtype) :: field_infra + + field_infra = handle%fields%get(field%label) + call read_field_chksum(field_infra, chksum, valid_chksum) +end subroutine read_field_chksum_infra + +!> Get the native (fieldtype) fields of a MOM framework file +subroutine get_file_fieldtypes(handle, fields) + class(MOM_infra_file), intent(in) :: handle + type(fieldtype), intent(out) :: fields(:) + + type(field_node_infra), pointer :: node + integer :: i + + ! NOTE: The tail is a pre-allocated empty node, so we check node%next + node => handle%fields%head + do i = 1, size(fields) + if (.not. associated(node%next)) & + call MOM_error(FATAL, 'fields(:) size exceeds number of registered fields.') + fields(i) = node%field + node => node%next + enddo +end subroutine get_file_fieldtypes + + +! MOM_netcdf_file methods + +!> Open a MOM netCDF file +subroutine open_file_nc(handle, filename, action, MOM_domain, threading, fileset) + class(MOM_netcdf_file), intent(inout) :: handle + character(len=*), intent(in) :: filename + integer, intent(in), optional :: action + type(MOM_domain_type), optional, intent(in) :: MOM_domain + integer, intent(in), optional :: threading + integer, intent(in), optional :: fileset + + if (.not. is_root_PE()) return + + call open_netcdf_file(handle%handle_nc, filename, action) + + handle%is_open = .true. + call handle%axes%init() + call handle%fields%init() +end subroutine open_file_nc + + +!> Close a MOM netCDF file +subroutine close_file_nc(handle) + class(MOM_netcdf_file), intent(inout) :: handle + + if (.not. is_root_PE()) return + + handle%is_open = .false. + call close_netcdf_file(handle%handle_nc) +end subroutine close_file_nc + + +!> Flush the buffer of a MOM netCDF file +subroutine flush_file_nc(handle) + class(MOM_netcdf_file), intent(in) :: handle + + if (.not. is_root_PE()) return + + call flush_netcdf_file(handle%handle_nc) +end subroutine flush_file_nc + + +!> Register an axis to the MOM netcdf file +function register_axis_nc(handle, label, units, longname, cartesian, sense, & + domain, data, edge_axis, calendar) result(axis) + class(MOM_netcdf_file), intent(inout) :: handle + !< Handle for a netCDF file that is open for writing + character(len=*), intent(in) :: label + !< The name in the file of this axis + character(len=*), intent(in) :: units + !< The units of this axis + character(len=*), intent(in) :: longname + !< The long description of this axis + character(len=*), optional, intent(in) :: cartesian + !< A variable indicating which direction this axis corresponds with. + !! Valid values include 'X', 'Y', 'Z', 'T', and 'N' for none. + integer, optional, intent(in) :: sense + !< This is 1 for axes whose values increase upward, or -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 + type(MOM_axis) :: axis + + type(netcdf_axis) :: axis_nc + + if (is_root_PE()) then + axis_nc = register_netcdf_axis(handle%handle_nc, label, units, longname, & + data, cartesian, sense) + + call handle%axes%append(axis_nc, label) + endif + axis%label = label +end function register_axis_nc + + +!> Register a field to the MOM netcdf file +function register_field_nc(handle, axes, label, units, longname, pack, & + standard_name, checksum) result(field) + class(MOM_netcdf_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_axis), intent(in) :: axes(:) + !< Handles for the axis used for this variable + character(len=*), intent(in) :: label + !< 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 + integer, optional, intent(in) :: pack + !< A precision reduction factor with which the variable. The default, 1, + !! has no reduction, but 2 is not uncommon. + character(len=*), optional, intent(in) :: standard_name + !< The standard (e.g., CMOR) name for this variable + integer(kind=int64), dimension(:), optional, intent(in) :: checksum + !< Checksum values that can be used to verify reads. + type(MOM_field) :: field + + type(netcdf_field) :: field_nc + type(netcdf_axis), allocatable :: axes_nc(:) + integer :: i + + if (is_root_PE()) then + allocate(axes_nc(size(axes))) + do i = 1, size(axes) + axes_nc(i) = handle%axes%get(axes(i)%label) + enddo + + field_nc = register_netcdf_field(handle%handle_nc, label, axes_nc, longname, units) + + call handle%fields%append(field_nc, label) + endif + field%label = label +end function register_field_nc + + +!> Write global metadata to the MOM netcdf file +subroutine write_attribute_nc(handle, name, attribute) + class(MOM_netcdf_file), intent(in) :: handle + !< Handle for a file that is open for writing + character(len=*), intent(in) :: name + !< The name in the file of this global attribute + character(len=*), intent(in) :: attribute + !< The value of this attribute + + if (.not. is_root_PE()) return + + call write_netcdf_attribute(handle%handle_nc, name, attribute) +end subroutine write_attribute_nc + + +!> Write a 4D field to the MOM netcdf file +subroutine write_field_4d_nc(handle, field_md, MOM_domain, field, tstamp, & + tile_count, fill_value) + class(MOM_netcdf_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_field), intent(in) :: field_md + !< Field type with metadata + type(MOM_domain_type), intent(in) :: MOM_domain + !< The MOM_Domain that describes the decomposition + real, intent(inout) :: field(:,:,:,:) + !< Field to write + 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 + + type(netcdf_field) :: field_nc + + if (.not. is_root_PE()) return + + field_nc = handle%fields%get(field_md%label) + call write_netcdf_field(handle%handle_nc, field_nc, field, time=tstamp) +end subroutine write_field_4d_nc + + +!> Write a 3D field to the MOM netcdf file +subroutine write_field_3d_nc(handle, field_md, MOM_domain, field, tstamp, & + tile_count, fill_value) + class(MOM_netcdf_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_field), intent(in) :: field_md + !< Field type with metadata + type(MOM_domain_type), intent(in) :: MOM_domain + !< The MOM_Domain that describes the decomposition + real, intent(inout) :: field(:,:,:) + !< Field to write + 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 + + type(netcdf_field) :: field_nc + + if (.not. is_root_PE()) return + + field_nc = handle%fields%get(field_md%label) + call write_netcdf_field(handle%handle_nc, field_nc, field, time=tstamp) +end subroutine write_field_3d_nc + + +!> Write a 2D field to the MOM netcdf file +subroutine write_field_2d_nc(handle, field_md, MOM_domain, field, tstamp, & + tile_count, fill_value) + class(MOM_netcdf_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_field), 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 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 + + type(netcdf_field) :: field_nc + + if (.not. is_root_PE()) return + + field_nc = handle%fields%get(field_md%label) + call write_netcdf_field(handle%handle_nc, field_nc, field, time=tstamp) +end subroutine write_field_2d_nc + + +!> Write a 1D field to the MOM netcdf file +subroutine write_field_1d_nc(handle, field_md, field, tstamp) + class(MOM_netcdf_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_field), 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 + + type(netcdf_field) :: field_nc + + if (.not. is_root_PE()) return + + field_nc = handle%fields%get(field_md%label) + call write_netcdf_field(handle%handle_nc, field_nc, field, time=tstamp) +end subroutine write_field_1d_nc + + +!> Write a 0D field to the MOM netcdf file +subroutine write_field_0d_nc(handle, field_md, field, tstamp) + class(MOM_netcdf_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_field), 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 + + type(netcdf_field) :: field_nc + + if (.not. is_root_PE()) return + + field_nc = handle%fields%get(field_md%label) + call write_netcdf_field(handle%handle_nc, field_nc, field, time=tstamp) +end subroutine write_field_0d_nc + + +!> Write an axis field to the MOM netcdf file +subroutine write_field_axis_nc(handle, axis) + class(MOM_netcdf_file), intent(inout) :: handle + !< Handle for a file that is open for writing + type(MOM_axis), intent(in) :: axis + !< An axis type variable with information to write + + type(netcdf_axis) :: axis_nc + + if (.not. is_root_PE()) return + + axis_nc = handle%axes%get(axis%label) + call write_netcdf_axis(handle%handle_nc, axis_nc) +end subroutine write_field_axis_nc + + +!> True if the framework file has been opened +logical function file_is_open_nc(handle) + class(MOM_netcdf_file), intent(in) :: handle + !< Handle to a file to inquire about + + file_is_open_nc = handle%is_open +end function file_is_open_nc + + +!> Return number of dimensions, variables, or time levels in a MOM netcdf file +subroutine get_file_info_nc(handle, ndim, nvar, ntime) + class(MOM_netcdf_file), intent(in) :: 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) :: ntime + !< The number of time levels in the file + + integer :: ndim_nc, nvar_nc + + if (.not. is_root_PE()) return + + call get_netcdf_size(handle%handle_nc, ndims=ndim_nc, nvars=nvar_nc, nsteps=ntime) + + ! MOM I/O follows legacy FMS behavior and excludes axes from field count + if (present(ndim)) ndim = ndim_nc + if (present(nvar)) nvar = nvar_nc - ndim_nc +end subroutine get_file_info_nc + + +!> Return the field metadata associated with a MOM netCDF file +subroutine get_file_fields_nc(handle, fields) + class(MOM_netcdf_file), intent(inout) :: handle + !< Handle for a file that is open for I/O + type(MOM_field), intent(inout) :: fields(:) + !< Field-type descriptions of all of the variables in a file. + + type(netcdf_axis), allocatable :: axes_nc(:) + type(netcdf_field), allocatable :: fields_nc(:) + integer :: i + + if (.not. is_root_PE()) return + + call get_netcdf_fields(handle%handle_nc, axes_nc, fields_nc) + if (size(fields) /= size(fields_nc)) & + call MOM_error(FATAL, 'Number of fields in file does not match field(:).') + + do i = 1, size(axes_nc) + call handle%axes%append(axes_nc(i), axes_nc(i)%label) + enddo + + do i = 1, size(fields) + fields(i)%label = trim(fields_nc(i)%label) + call handle%fields%append(fields_nc(i), fields_nc(i)%label) + enddo +end subroutine get_file_fields_nc + + +!> Get attributes from a netCDF field +subroutine get_field_atts_nc(handle, field, name, units, longname, checksum) + class(MOM_netcdf_file), intent(in) :: handle + !< File where field is stored + type(MOM_field), intent(in) :: field + !< The field to extract information from + character(len=*), optional, intent(out) :: name + !< The variable name + character(len=*), optional, intent(out) :: units + !< The units of the variable + character(len=*), optional, intent(out) :: longname + !< The long name of the variable + integer(kind=int64), optional, intent(out) :: checksum(:) + !< The checksums of the variable in a file + + call MOM_error(FATAL, 'get_field_atts over netCDF is not yet implemented.') +end subroutine get_field_atts_nc + + +!> Interface to read_field_chksum +subroutine read_field_chksum_nc(handle, field, chksum, valid_chksum) + class(MOM_netcdf_file), intent(in) :: handle + !< File where field is stored + type(MOM_field), 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 + + call MOM_error(FATAL, 'read_field_chksum over netCDF is not yet implemented.') +end subroutine read_field_chksum_nc + + +!> \namespace MOM_IO_file +!! +!! This file defines the MOM_file classes used to inferface with the internal +!! IO handlers, such as the configured "infra" layer (FMS) or native netCDF. +!! +!! `MOM_file`: The generic class used to reference any file type +!! Cannot be used in a variable declaration. +!! +!! `MOM_infra_file`: A file handler for use by the infra layer. Currently this +!! means an FMS file, such a restart or diagnostic output. +!! +!! `MOM_netcdf_file`: A netCDF file handler for MOM-specific I/O. This may +!! include operations outside the scope of FMS or other infra frameworks. + +end module MOM_io_file diff --git a/src/framework/MOM_netcdf.F90 b/src/framework/MOM_netcdf.F90 new file mode 100644 index 0000000000..d09ae5cf95 --- /dev/null +++ b/src/framework/MOM_netcdf.F90 @@ -0,0 +1,755 @@ +!> MOM6 interface to netCDF operations +module MOM_netcdf + +! This file is part of MOM6. See LICENSE.md for the license. + +use, intrinsic :: iso_fortran_env, only : real32, real64 + +use netcdf, only : nf90_create, nf90_open, nf90_close +use netcdf, only : nf90_sync +use netcdf, only : NF90_CLOBBER, NF90_NOCLOBBER, NF90_WRITE, NF90_NOWRITE +use netcdf, only : nf90_enddef +use netcdf, only : nf90_def_dim, nf90_def_var +use netcdf, only : NF90_UNLIMITED +use netcdf, only : nf90_get_var +use netcdf, only : nf90_put_var, nf90_put_att +use netcdf, only : NF90_FLOAT, NF90_DOUBLE +use netcdf, only : nf90_strerror, NF90_NOERR +use netcdf, only : NF90_GLOBAL +use netcdf, only : nf90_inquire, nf90_inquire_dimension, nf90_inquire_variable +use netcdf, only : nf90_inq_dimids, nf90_inq_varids +use netcdf, only : NF90_MAX_NAME + +use MOM_error_handler, only : MOM_error, FATAL +use MOM_io_infra, only : READONLY_FILE, WRITEONLY_FILE +use MOM_io_infra, only : APPEND_FILE, OVERWRITE_FILE + +implicit none ; private + +public :: netcdf_file_type +public :: netcdf_axis +public :: netcdf_field +public :: open_netcdf_file +public :: close_netcdf_file +public :: flush_netcdf_file +public :: register_netcdf_axis +public :: register_netcdf_field +public :: write_netcdf_field +public :: write_netcdf_axis +public :: write_netcdf_attribute +public :: get_netcdf_size +public :: get_netcdf_fields + + +!> Internal time value used to indicate an uninitialized time +real, parameter :: NULLTIME = -1 +! NOTE: For now, we use the FMS-compatible value, but may change in the future. + + +!> netCDF file abstraction +type :: netcdf_file_type + private + integer :: ncid + !< netCDF file ID + character(len=:), allocatable :: filename + !< netCDF filename + logical :: define_mode + !< True if file is in define mode. + integer :: time_id + !< Time axis variable ID + real :: time + !< Current model time + integer :: time_level + !< Current time level for output +end type netcdf_file_type + + +!> Dimension axis for a netCDF file +type :: netcdf_axis + private + character(len=:), allocatable, public :: label + !< Axis label name + real, allocatable :: points(:) + !< Grid points along the axis + integer :: dimid + !< netCDF dimension ID associated with axis + integer :: varid + !< netCDF variable ID associated with axis +end type netcdf_axis + + +!> Field variable for a netCDF file +type netcdf_field + private + character(len=:), allocatable, public :: label + !< Variable name + integer :: varid + !< netCDF variable ID for field +end type netcdf_field + + +!> Write values to a field of a netCDF file +interface write_netcdf_field + module procedure write_netcdf_field_4d + module procedure write_netcdf_field_3d + module procedure write_netcdf_field_2d + module procedure write_netcdf_field_1d + module procedure write_netcdf_field_0d +end interface write_netcdf_field + +contains + +subroutine open_netcdf_file(handle, filename, mode) + type(netcdf_file_type), intent(inout) :: handle + !< netCDF file handle + character(len=*), intent(in) :: filename + !< netCDF filename + integer, intent(in), optional :: mode + !< Input MOM I/O mode + + integer :: io_mode + ! MOM I/O mode + integer :: cmode + ! netCDF creation mode + integer :: rc + ! nf90_create return code + character(len=:), allocatable :: msg + ! netCDF error message buffer + + ! I/O configuration + io_mode = WRITEONLY_FILE + if (present(mode)) io_mode = mode + + ! Translate the MOM I/O config to the netCDF mode + select case(io_mode) + case (WRITEONLY_FILE) + rc = nf90_create(filename, nf90_noclobber, handle%ncid) + handle%define_mode = .true. + case (OVERWRITE_FILE) + rc = nf90_create(filename, nf90_clobber, handle%ncid) + handle%define_mode = .true. + case (APPEND_FILE) + rc = nf90_open(filename, nf90_write, handle%ncid) + handle%define_mode = .false. + case (READONLY_FILE) + rc = nf90_open(filename, nf90_nowrite, handle%ncid) + handle%define_mode = .false. + case default + call MOM_error(FATAL, & + 'open_netcdf_file: File ' // filename // ': Unknown mode.') + end select + call check_netcdf_call(rc, 'open_netcdf_file', 'File ' // filename) + + handle%filename = filename + + ! FMS writes the filename as an attribute + call write_netcdf_attribute(handle, 'filename', filename) +end subroutine open_netcdf_file + + +!> Close an opened netCDF file. +subroutine close_netcdf_file(handle) + type(netcdf_file_type), intent(in) :: handle + + integer :: rc + + rc = nf90_close(handle%ncid) + call check_netcdf_call(rc, 'close_netcdf_file', & + 'File "' // handle%filename // '"') +end subroutine close_netcdf_file + + +!> Flush buffered output to the netCDF file +subroutine flush_netcdf_file(handle) + type(netcdf_file_type), intent(in) :: handle + + integer :: rc + + rc = nf90_sync(handle%ncid) + call check_netcdf_call(rc, 'flush_netcdf_file', & + 'File "' // handle%filename // '"') +end subroutine flush_netcdf_file + + +!> Change netCDF mode of handle from 'define' to 'write'. +subroutine enable_netcdf_write(handle) + type(netcdf_file_type), intent(inout) :: handle + + integer :: rc + + if (handle%define_mode) then + rc = nf90_enddef(handle%ncid) + call check_netcdf_call(rc, 'enable_netcdf_write', & + 'File "' // handle%filename // '"') + handle%define_mode = .false. + endif +end subroutine enable_netcdf_write + + +!> Register a netCDF variable +function register_netcdf_field(handle, label, axes, longname, units) & + result(field) + type(netcdf_file_type), intent(in) :: handle + !< netCDF file handle + character(len=*), intent(in) :: label + !< netCDF field name in the file + type(netcdf_axis), intent(in) :: axes(:) + !< Axes along which field is defined + character(len=*), intent(in) :: longname + !< Long name of the netCDF field + character(len=*), intent(in) :: units + !< Field units of measurement + type(netcdf_field) :: field + !< netCDF field + + integer :: rc + ! netCDF function return code + integer :: i + ! Loop index + integer, allocatable :: dimids(:) + ! netCDF dimension IDs of axes + integer :: xtype + ! netCDF data type + + ! Gather the axis netCDF dimension IDs + allocate(dimids(size(axes))) + dimids(:) = [(axes(i)%dimid, i = 1, size(axes))] + + ! Determine the corresponding netCDF data type + ! TODO: Support a `pack`-like argument + select case (kind(1.0)) + case (real32) + xtype = NF90_FLOAT + case (real64) + xtype = NF90_DOUBLE + case default + call MOM_error(FATAL, "register_netcdf_axis: Unknown kind(real).") + end select + + ! Register the field variable + rc = nf90_def_var(handle%ncid, label, xtype, dimids, field%varid) + call check_netcdf_call(rc, 'register_netcdf_field', & + 'File "' // handle%filename // '", Field "' // label // '"') + + ! Assign attributes + + rc = nf90_put_att(handle%ncid, field%varid, 'long_name', longname) + call check_netcdf_call(rc, 'register_netcdf_field', & + 'Attribute "long_name" of variable "' // label // '" in file "' & + // handle%filename // '"') + + rc = nf90_put_att(handle%ncid, field%varid, 'units', units) + call check_netcdf_call(rc, 'register_netcdf_field', & + 'Attribute "units" of variable "' // label // '" in file "' & + // handle%filename // '"') +end function register_netcdf_field + + +!> Create an axis and associated dimension in a netCDF file +function register_netcdf_axis(handle, label, units, longname, points, & + cartesian, sense) result(axis) + type(netcdf_file_type), intent(inout) :: handle + !< netCDF file handle + character(len=*), intent(in) :: label + !< netCDF axis name in the file + character(len=*), intent(in), optional :: units + !< Axis units of measurement + character(len=*), intent(in), optional :: longname + !< Long name of the axis + real, intent(in), optional :: points(:) + !< Values of axis points (for fixed axes) + character(len=*), intent(in), optional :: cartesian + !< Character denoting axis direction: X, Y, Z, T, or N for none + integer, intent(in), optional :: sense + !< Axis direction; +1 if axis increases upward or -1 if downward + + type(netcdf_axis) :: axis + !< netCDF coordinate axis + + integer :: xtype + ! netCDF external data type + integer :: rc + ! netCDF function return code + logical :: unlimited + ! True if the axis is unlimited in size (e.g. time) + integer :: axis_size + ! Either the number of points in the axis, or unlimited flag + integer :: axis_sense + ! Axis direction; +1 if axis increases upward or -1 if downward + character(len=:), allocatable :: sense_attr + ! CF-compiant value of sense attribute (as 'positive') + + ! Create the axis dimension + unlimited = .false. + if (present(cartesian)) then + if (cartesian == 'T') unlimited = .true. + endif + + ! Either the axis is explicitly set with data or is declared as unlimited + if (present(points) .eqv. unlimited) then + call MOM_error(FATAL, & + "Axis must either have explicit points or be a time axis ('T').") + endif + + if (present(points)) then + axis_size = size(points) + allocate(axis%points(axis_size)) + axis%points(:) = points(:) + else + axis_size = NF90_UNLIMITED + endif + + rc = nf90_def_dim(handle%ncid, label, axis_size, axis%dimid) + call check_netcdf_call(rc, 'register_netcdf_axis', & + 'Dimension "' // label // '" in file "' // handle%filename // '"') + + ! Determine the corresponding netCDF data type + ! TODO: Support a `pack`-like argument + select case (kind(1.0)) + case (real32) + xtype = NF90_FLOAT + case (real64) + xtype = NF90_DOUBLE + case default + call MOM_error(FATAL, "register_netcdf_axis: Unknown kind(real).") + end select + + ! Create a variable corresponding to the axis + rc = nf90_def_var(handle%ncid, label, xtype, axis%dimid, axis%varid) + call check_netcdf_call(rc, 'register_netcdf_axis', & + 'Variable ' // label // ' in file ' // handle%filename) + + ! Define the time axis, if available + if (unlimited) then + handle%time_id = axis%varid + handle%time_level = 0 + handle%time = NULLTIME + endif + + ! Assign attributes if present + if (present(longname)) then + rc = nf90_put_att(handle%ncid, axis%varid, 'long_name', longname) + call check_netcdf_call(rc, 'register_netcdf_axis', & + 'Attribute ''long_name'' of variable ' // label // ' in file ' & + // handle%filename) + endif + + if (present(units)) then + rc = nf90_put_att(handle%ncid, axis%varid, 'units', units) + call check_netcdf_call(rc, 'register_netcdf_axis', & + 'Attribute ''units'' of variable ' // label // ' in file ' & + // handle%filename) + endif + + if (present(cartesian)) then + rc = nf90_put_att(handle%ncid, axis%varid, 'cartesian_axis', cartesian) + call check_netcdf_call(rc, 'register_netcdf_axis', & + 'Attribute ''cartesian_axis'' of variable ' // label // ' in file ' & + // handle%filename) + endif + + axis_sense = 0 + if (present(sense)) axis_sense = sense + + if (axis_sense /= 0) then + select case (axis_sense) + case (1) + sense_attr = 'up' + case (-1) + sense_attr = 'down' + case default + call MOM_error(FATAL, 'register_netcdf_axis: sense must be either ' & + // '0, 1, or -1.') + end select + rc = nf90_put_att(handle%ncid, axis%varid, 'positive', sense_attr) + call check_netcdf_call(rc, 'register_netcdf_axis', & + 'Attribute "positive" of variable "' // label // '" in file "' & + // handle%filename // '"') + endif +end function register_netcdf_axis + + +!> Write a 4D array to a compatible netCDF field +subroutine write_netcdf_field_4d(handle, field, values, time) + type(netcdf_file_type), intent(inout) :: handle + !< netCDF file handle + type(netcdf_field), intent(in) :: field + !< Field metadata + real, intent(in) :: values(:,:,:,:) + !< Field values + real, intent(in), optional :: time + !< Timestep index to write data + + integer :: rc + ! netCDF return code + integer :: start(5) + ! Start indices, if timestep is included + + ! Verify write mode + if (handle%define_mode) & + call enable_netcdf_write(handle) + + if (present(time)) then + call update_netcdf_timestep(handle, time) + start(:4) = 1 + start(5) = handle%time_level + rc = nf90_put_var(handle%ncid, field%varid, values, start) + else + rc = nf90_put_var(handle%ncid, field%varid, values) + endif + call check_netcdf_call(rc, 'write_netcdf_file', & + 'File "' // handle%filename // '", Field "' // field%label // '"') +end subroutine write_netcdf_field_4d + + +!> Write a 3D array to a compatible netCDF field +subroutine write_netcdf_field_3d(handle, field, values, time) + type(netcdf_file_type), intent(inout) :: handle + !< netCDF file handle + type(netcdf_field), intent(in) :: field + !< Field metadata + real, intent(in) :: values(:,:,:) + !< Field values + real, intent(in), optional :: time + !< Timestep index to write data + + integer :: rc + ! netCDF return code + integer :: start(4) + ! Start indices, if timestep is included + + ! Verify write mode + if (handle%define_mode) & + call enable_netcdf_write(handle) + + if (present(time)) then + call update_netcdf_timestep(handle, time) + start(:3) = 1 + start(4) = handle%time_level + rc = nf90_put_var(handle%ncid, field%varid, values, start) + else + rc = nf90_put_var(handle%ncid, field%varid, values) + endif + call check_netcdf_call(rc, 'write_netcdf_file', & + 'File "' // handle%filename // '", Field "' // field%label // '"') +end subroutine write_netcdf_field_3d + + +!> Write a 2D array to a compatible netCDF field +subroutine write_netcdf_field_2d(handle, field, values, time) + type(netcdf_file_type), intent(inout) :: handle + !< netCDF file handle + type(netcdf_field), intent(in) :: field + !< Field metadata + real, intent(in) :: values(:,:) + !< Field values + real, intent(in), optional :: time + !< Timestep index to write data + + integer :: rc + ! netCDF return code + integer :: start(3) + ! Start indices, if timestep is included + + ! Verify write mode + if (handle%define_mode) & + call enable_netcdf_write(handle) + + if (present(time)) then + call update_netcdf_timestep(handle, time) + start(:2) = 1 + start(3) = handle%time_level + rc = nf90_put_var(handle%ncid, field%varid, values, start) + else + rc = nf90_put_var(handle%ncid, field%varid, values) + endif + call check_netcdf_call(rc, 'write_netcdf_file', & + 'File "' // handle%filename // '", Field "' // field%label // '"') +end subroutine write_netcdf_field_2d + + +!> Write a 1D array to a compatible netCDF field +subroutine write_netcdf_field_1d(handle, field, values, time) + type(netcdf_file_type), intent(inout) :: handle + !< netCDF file handle + type(netcdf_field), intent(in) :: field + !< Field metadata + real, intent(in) :: values(:) + !< Field values + real, intent(in), optional :: time + !< Timestep index to write data + + integer :: rc + ! netCDF return code + integer :: start(2) + ! Start indices, if timestep is included + + ! Verify write mode + if (handle%define_mode) & + call enable_netcdf_write(handle) + + if (present(time)) then + call update_netcdf_timestep(handle, time) + start(1) = 1 + start(2) = handle%time_level + rc = nf90_put_var(handle%ncid, field%varid, values, start) + else + rc = nf90_put_var(handle%ncid, field%varid, values) + endif + call check_netcdf_call(rc, 'write_netcdf_file', & + 'File "' // handle%filename // '", Field "' // field%label // '"') +end subroutine write_netcdf_field_1d + + +!> Write a scalar to a compatible netCDF field +subroutine write_netcdf_field_0d(handle, field, scalar, time) + type(netcdf_file_type), intent(inout) :: handle + !< netCDF file handle + type(netcdf_field), intent(in) :: field + !< Field metadata + real, intent(in) :: scalar + !< Field values + real, intent(in), optional :: time + !< Timestep index to write data + + integer :: rc + ! netCDF return code + integer :: start(1) + ! Start indices, if timestep is included + + ! Verify write mode + if (handle%define_mode) & + call enable_netcdf_write(handle) + + if (present(time)) then + call update_netcdf_timestep(handle, time) + start(1) = handle%time_level + rc = nf90_put_var(handle%ncid, field%varid, scalar, start) + else + rc = nf90_put_var(handle%ncid, field%varid, scalar) + endif + call check_netcdf_call(rc, 'write_netcdf_file', & + 'File "' // handle%filename // '", Field "' // field%label // '"') +end subroutine write_netcdf_field_0d + + +!> Write axis points to associated netCDF variable +subroutine write_netcdf_axis(handle, axis) + type(netcdf_file_type), intent(inout) :: handle + !< netCDF file handle + type(netcdf_axis), intent(in) :: axis + !< field variable + + integer :: rc + ! netCDF return code + + ! Verify write mode + if (handle%define_mode) & + call enable_netcdf_write(handle) + + rc = nf90_put_var(handle%ncid, axis%varid, axis%points) + call check_netcdf_call(rc, 'write_netcdf_axis', & + 'File "' // handle%filename // '", Axis "' // axis%label // '"') +end subroutine write_netcdf_axis + + +!> Write a global attribute to a netCDF file +subroutine write_netcdf_attribute(handle, label, attribute) + type(netcdf_file_type), intent(in) :: handle + !< netCDF file handle + character(len=*), intent(in) :: label + !< File attribute + character(len=*), intent(in) :: attribute + !< File attribute value + + integer :: rc + ! netCDF return code + + rc = nf90_put_att(handle%ncid, NF90_GLOBAL, label, attribute) + call check_netcdf_call(rc, 'write_netcdf_attribute', & + 'File "' // handle%filename // '", Attribute "' // label // '"') +end subroutine write_netcdf_attribute + + +! This is a thin interface to nf90_inquire, designed to mirror the existing +! I/O API. A more axis-aware system might not need this, but for now it's here +!> Get the number of dimensions, variables, and timesteps in a netCDF file +subroutine get_netcdf_size(handle, ndims, nvars, nsteps) + type(netcdf_file_type), intent(in) :: handle + !< netCDF input file + integer, intent(out), optional :: ndims + !< number of dimensions in the file + integer, intent(out), optional :: nvars + !< number of variables in the file + integer, intent(out), optional :: nsteps + !< number of values in the file's unlimited axis + + integer :: rc + ! netCDF return code + integer :: unlimited_dimid + ! netCDF dimension ID for unlimited time axis + + rc = nf90_inquire(handle%ncid, & + nDimensions=ndims, & + nVariables=nvars, & + unlimitedDimId=unlimited_dimid & + ) + call check_netcdf_call(rc, 'get_netcdf_size', & + 'File "' // handle%filename // '"') + + rc = nf90_inquire_dimension(handle%ncid, unlimited_dimid, len=nsteps) + call check_netcdf_call(rc, 'get_netcdf_size', & + 'File "' // handle%filename // '"') +end subroutine get_netcdf_size + + +!> Get the metadata of the registered fields in a netCDF file +subroutine get_netcdf_fields(handle, axes, fields) + type(netcdf_file_type), intent(inout) :: handle + type(netcdf_axis), intent(inout), allocatable :: axes(:) + type(netcdf_field), intent(inout), allocatable :: fields(:) + + integer :: ndims + ! Number of netCDF dimensions + integer :: nvars + ! Number of netCDF dimensions + integer :: nfields + ! Number of fields in the file (i.e. non-axis variables) + integer, allocatable :: dimids(:) + ! netCDF dimension IDs of file + integer, allocatable :: varids(:) + ! netCDF variable IDs of file + integer :: unlim_dimid + ! netCDF dimension ID for the unlimited axis variable, if present + integer :: unlim_index + ! Index of the unlimited axis in axes(:), if present + character(len=NF90_MAX_NAME) :: label + ! Current dimension or variable label + integer :: len + ! Current dimension length + integer :: rc + ! netCDF return code + integer :: grp_ndims, grp_nvars + ! Group-based counts for nf90_inq_* (unused) + logical :: is_axis + ! True if the current variable is an axis + integer :: i, j, n + + integer, save :: no_parent_groups = 0 + ! Flag indicating exclusion of parent groups in netCDF file + ! NOTE: This must be passed as a variable, and cannot be declared as a + ! parameter. + + rc = nf90_inquire(handle%ncid, & + nDimensions=ndims, & + nVariables=nvars, & + unlimitedDimId=unlim_dimid & + ) + call check_netcdf_call(rc, 'get_netcdf_fields', & + 'File "' // handle%filename // '"') + + allocate(dimids(ndims)) + rc = nf90_inq_dimids(handle%ncid, grp_ndims, dimids, no_parent_groups) + call check_netcdf_call(rc, 'get_netcdf_fields', & + 'File "' // handle%filename // '"') + + allocate(varids(nvars)) + rc = nf90_inq_varids(handle%ncid, grp_nvars, varids) + call check_netcdf_call(rc, 'get_netcdf_fields', & + 'File "' // handle%filename // '"') + + allocate(axes(ndims)) + do i = 1, ndims + rc = nf90_inquire_dimension(handle%ncid, dimids(i), name=label, len=len) + call check_netcdf_call(rc, 'get_netcdf_fields', & + 'File "' // handle%filename // '"') + + ! Check for the unlimited axis + if (dimids(i) == unlim_dimid) unlim_index = i + + axes(i)%dimid = dimids(i) + axes(i)%label = trim(label) + allocate(axes(i)%points(len)) + enddo + + nfields = nvars - ndims + allocate(fields(nfields)) + + n = 0 + do i = 1, nvars + rc = nf90_inquire_variable(handle%ncid, varids(i), name=label) + call check_netcdf_call(rc, 'get_netcdf_fields', & + 'File "' // handle%filename // '"') + + ! Check if variable is an axis + is_axis = .false. + do j = 1, ndims + if (label == axes(j)%label) then + rc = nf90_get_var(handle%ncid, varids(i), axes(j)%points) + call check_netcdf_call(rc, 'get_netcdf_fields', & + 'File "' // handle%filename // '"') + axes(j)%varid = varids(i) + + if (j == unlim_index) then + handle%time_id = varids(i) + handle%time_level = size(axes(j)%points) + handle%time = NULLTIME + endif + + is_axis = .true. + exit + endif + enddo + if (is_axis) cycle + + n = n + 1 + fields(n)%label = trim(label) + fields(n)%varid = varids(i) + enddo +end subroutine get_netcdf_fields + + +subroutine update_netcdf_timestep(handle, time) + type(netcdf_file_type), intent(inout) :: handle + !< netCDF file handle + real, intent(in) :: time + !< New model time + + integer :: start(1) + !< Time axis start index array + integer :: rc + !< netCDF return code + + if (time > handle%time + epsilon(time)) then + handle%time = time + handle%time_level = handle%time_level + 1 + + ! Write new value to time axis + start = [handle%time_level] + rc = nf90_put_var(handle%ncid, handle%time_id, time, start=start) + call check_netcdf_call(rc, 'update_netcdf_timestep', & + 'File "' // handle%filename // '"') + endif +end subroutine update_netcdf_timestep + + +!> Check netCDF function return codes, report the error log, and abort the run. +subroutine check_netcdf_call(ncerr, header, message) + integer, intent(in) :: ncerr + !< netCDF error code + character(len=*), intent(in) :: header + !< Message header (usually calling subroutine) + character(len=*), intent(in) :: message + !< Error message (usually action which instigated the error) + + character(len=:), allocatable :: errmsg + ! Full error message, including netCDF message + + if (ncerr /= nf90_noerr) then + errmsg = trim(header) // ": " // trim(message) // new_line('/') & + // trim(nf90_strerror(ncerr)) + call MOM_error(FATAL, errmsg) + endif +end subroutine check_netcdf_call + +end module MOM_netcdf diff --git a/src/framework/MOM_restart.F90 b/src/framework/MOM_restart.F90 index a76e96499f..2939a7d907 100644 --- a/src/framework/MOM_restart.F90 +++ b/src/framework/MOM_restart.F90 @@ -8,9 +8,9 @@ module MOM_restart use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, is_root_pe use MOM_file_parser, only : get_param, log_param, log_version, param_file_type use MOM_grid, only : ocean_grid_type -use MOM_io, only : create_file, file_type, fieldtype, file_exists, open_file, close_file -use MOM_io, only : MOM_read_data, read_data, MOM_write_field, read_field_chksum, field_exists -use MOM_io, only : get_file_info, get_file_fields, get_field_atts, get_file_times +use MOM_io, only : create_MOM_file, file_exists +use MOM_io, only : MOM_infra_file, MOM_field +use MOM_io, only : MOM_read_data, read_data, MOM_write_field, field_exists use MOM_io, only : vardesc, var_desc, query_vardesc, modify_vardesc, get_filename_appendix use MOM_io, only : MULTIPLE, READONLY_FILE, SINGLE_FILE use MOM_io, only : CENTER, CORNER, NORTH_FACE, EAST_FACE @@ -1258,7 +1258,7 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ ! Local variables type(vardesc) :: vars(CS%max_fields) ! Descriptions of the fields that ! are to be read from the restart file. - type(fieldtype) :: fields(CS%max_fields) ! Opaque types containing metadata describing + type(MOM_field) :: fields(CS%max_fields) ! Opaque types containing metadata describing ! each variable that will be written. character(len=512) :: restartpath ! The restart file path (dir/file). character(len=256) :: restartname ! The restart file name (no dir). @@ -1272,7 +1272,7 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ ! versions of NetCDF, the value was 2147483647_8. integer :: start_var, next_var ! The starting variables of the ! current and next files. - type(file_type) :: IO_handle ! The I/O handle of the open fileset + type(MOM_infra_file) :: IO_handle ! The I/O handle of the open fileset integer :: m, nz integer :: num_files ! The number of restart files that will be used. integer :: seconds, days, year, month, hour, minute @@ -1408,11 +1408,11 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ enddo if (CS%parallel_restartfiles) then - call create_file(IO_handle, trim(restartpath), vars, (next_var-start_var), & - fields, MULTIPLE, G=G, GV=GV, checksums=check_val) + call create_MOM_file(IO_handle, trim(restartpath), vars, next_var-start_var, & + fields, MULTIPLE, G=G, GV=GV, checksums=check_val) else - call create_file(IO_handle, trim(restartpath), vars, (next_var-start_var), & - fields, SINGLE_FILE, G=G, GV=GV, checksums=check_val) + call create_MOM_file(IO_handle, trim(restartpath), vars, next_var-start_var, & + fields, SINGLE_FILE, G=G, GV=GV, checksums=check_val) endif do m=start_var,next_var-1 @@ -1434,7 +1434,7 @@ subroutine save_restart(directory, time, G, CS, time_stamped, filename, GV, num_ endif enddo - call close_file(IO_handle) + call IO_handle%close() num_files = num_files+1 @@ -1466,14 +1466,14 @@ subroutine restore_state(filename, directory, day, G, CS) integer :: isL, ieL, jsL, jeL integer :: nvar, ntime, pos - type(file_type) :: IO_handles(CS%max_fields) ! The I/O units of all open files. + type(MOM_infra_file) :: IO_handles(CS%max_fields) ! The I/O units of all open files. character(len=200) :: unit_path(CS%max_fields) ! The file names. logical :: unit_is_global(CS%max_fields) ! True if the file is global. character(len=8) :: hor_grid ! Variable grid info. real :: t1, t2 ! Two times. real, allocatable :: time_vals(:) - type(fieldtype), allocatable :: fields(:) + type(MOM_field), allocatable :: fields(:) logical :: is_there_a_checksum ! Is there a valid checksum that should be checked. integer(kind=8) :: checksum_file ! The checksum value recorded in the input file. integer(kind=8) :: checksum_data ! The checksum value for the data that was read in. @@ -1500,7 +1500,7 @@ subroutine restore_state(filename, directory, day, G, CS) ! Get the time from the first file in the list that has one. do n=1,num_file - call get_file_times(IO_handles(n), time_vals, ntime) + call IO_handles(n)%get_file_times(time_vals, ntime) if (ntime < 1) cycle t1 = time_vals(1) @@ -1516,7 +1516,7 @@ subroutine restore_state(filename, directory, day, G, CS) ! Check the remaining files for different times and issue a warning ! if they differ from the first time. do m = n+1,num_file - call get_file_times(IO_handles(n), time_vals, ntime) + call IO_handles(n)%get_file_times(time_vals, ntime) if (ntime < 1) cycle t2 = time_vals(1) @@ -1532,13 +1532,13 @@ subroutine restore_state(filename, directory, day, G, CS) ! Read each variable from the first file in which it is found. do n=1,num_file - call get_file_info(IO_handles(n), nvar=nvar) + call IO_handles(n)%get_file_info(nvar=nvar) allocate(fields(nvar)) - call get_file_fields(IO_handles(n), fields(1:nvar)) + call IO_handles(n)%get_file_fields(fields(1:nvar)) do m=1, nvar - call get_field_atts(fields(m), name=varname) + call IO_handles(n)%get_field_atts(fields(m), name=varname) do i=1,CS%num_obsolete_vars if (adjustl(lowercase(trim(varname))) == adjustl(lowercase(trim(CS%restart_obsolete(i)%field_name)))) then call MOM_error(FATAL, "MOM_restart restore_state: Attempting to use obsolete restart field "//& @@ -1571,11 +1571,11 @@ subroutine restore_state(filename, directory, day, G, CS) call get_checksum_loop_ranges(G, pos, isL, ieL, jsL, jeL) do i=1, nvar - call get_field_atts(fields(i), name=varname) + call IO_handles(n)%get_field_atts(fields(i), name=varname) if (lowercase(trim(varname)) == lowercase(trim(CS%restart_field(m)%var_name))) then checksum_data = -1 if (CS%checksum_required) then - call read_field_chksum(fields(i), checksum_file, is_there_a_checksum) + call IO_handles(n)%read_field_chksum(fields(i), checksum_file, is_there_a_checksum) else checksum_file = -1 is_there_a_checksum = .false. ! Do not need to do data checksumming. @@ -1643,7 +1643,7 @@ subroutine restore_state(filename, directory, day, G, CS) enddo do n=1,num_file - call close_file(IO_handles(n)) + call IO_handles(n)%close() enddo ! Check whether any mandatory fields have not been found. @@ -1745,7 +1745,7 @@ function open_restart_units(filename, directory, G, CS, IO_handles, file_paths, type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure type(MOM_restart_CS), intent(in) :: CS !< MOM restart control struct - type(file_type), dimension(:), & + type(MOM_infra_file), dimension(:), & optional, intent(out) :: IO_handles !< The I/O handles of all opened files character(len=*), dimension(:), & optional, intent(out) :: file_paths !< The full paths to open files @@ -1822,7 +1822,7 @@ function open_restart_units(filename, directory, G, CS, IO_handles, file_paths, if (fexists) then nf = nf + 1 if (present(IO_handles)) & - call open_file(IO_handles(nf), trim(filepath), READONLY_FILE, & + call IO_handles(nf)%open(trim(filepath), READONLY_FILE, & threading=MULTIPLE, fileset=SINGLE_FILE) if (present(global_files)) global_files(nf) = .true. if (present(file_paths)) file_paths(nf) = filepath @@ -1832,7 +1832,7 @@ function open_restart_units(filename, directory, G, CS, IO_handles, file_paths, if (fexists) then nf = nf + 1 if (present(IO_handles)) & - call open_file(IO_handles(nf), trim(filepath), READONLY_FILE, MOM_domain=G%Domain) + call IO_handles(nf)%open(trim(filepath), READONLY_FILE, MOM_domain=G%Domain) if (present(global_files)) global_files(nf) = .false. if (present(file_paths)) file_paths(nf) = filepath endif @@ -1854,7 +1854,7 @@ function open_restart_units(filename, directory, G, CS, IO_handles, file_paths, if (fexists) then nf = nf + 1 if (present(IO_handles)) & - call open_file(IO_handles(nf), trim(filepath), READONLY_FILE, & + call IO_handles(nf)%open(trim(filepath), READONLY_FILE, & threading=MULTIPLE, fileset=SINGLE_FILE) if (present(global_files)) global_files(nf) = .true. if (present(file_paths)) file_paths(nf) = filepath diff --git a/src/ice_shelf/MOM_ice_shelf.F90 b/src/ice_shelf/MOM_ice_shelf.F90 index aaa53dee59..a78c17803c 100644 --- a/src/ice_shelf/MOM_ice_shelf.F90 +++ b/src/ice_shelf/MOM_ice_shelf.F90 @@ -33,7 +33,7 @@ module MOM_ice_shelf use user_initialization, only : user_initialize_topography use MOM_io, only : field_exists, file_exists, MOM_read_data, write_version_number use MOM_io, only : slasher, fieldtype, vardesc, var_desc -use MOM_io, only : write_field, close_file, SINGLE_FILE, MULTIPLE +use MOM_io, only : close_file, SINGLE_FILE, MULTIPLE use MOM_restart, only : register_restart_field, save_restart use MOM_restart, only : restart_init, restore_state, MOM_restart_CS, register_restart_pair use MOM_time_manager, only : time_type, time_type_to_real, real_to_time, operator(>), operator(-) diff --git a/src/initialization/MOM_coord_initialization.F90 b/src/initialization/MOM_coord_initialization.F90 index c956c6b4be..78f739c461 100644 --- a/src/initialization/MOM_coord_initialization.F90 +++ b/src/initialization/MOM_coord_initialization.F90 @@ -8,7 +8,8 @@ module MOM_coord_initialization use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, read_param, log_param, param_file_type, log_version -use MOM_io, only : close_file, create_file, file_type, fieldtype, file_exists +use MOM_io, only : create_MOM_file, file_exists +use MOM_io, only : MOM_infra_file, MOM_field use MOM_io, only : MOM_read_data, MOM_write_field, vardesc, var_desc, SINGLE_FILE use MOM_string_functions, only : slasher, uppercase use MOM_unit_scaling, only : unit_scale_type @@ -526,20 +527,21 @@ subroutine write_vertgrid_file(GV, US, param_file, directory) ! Local variables character(len=240) :: filepath type(vardesc) :: vars(2) - type(fieldtype) :: fields(2) - type(file_type) :: IO_handle ! The I/O handle of the fileset + type(MOM_field) :: fields(2) + type(MOM_infra_file) :: IO_handle ! The I/O handle of the fileset filepath = trim(directory) // trim("Vertical_coordinate") vars(1) = var_desc("R","kilogram meter-3","Target Potential Density",'1','L','1') vars(2) = var_desc("g","meter second-2","Reduced gravity",'1','L','1') - call create_file(IO_handle, trim(filepath), vars, 2, fields, SINGLE_FILE, GV=GV) + call create_MOM_file(IO_handle, trim(filepath), vars, 2, fields, & + SINGLE_FILE, GV=GV) call MOM_write_field(IO_handle, fields(1), GV%Rlay, scale=US%R_to_kg_m3) call MOM_write_field(IO_handle, fields(2), GV%g_prime, scale=US%L_T_to_m_s**2*US%m_to_Z) - call close_file(IO_handle) + call IO_handle%close() end subroutine write_vertgrid_file diff --git a/src/initialization/MOM_shared_initialization.F90 b/src/initialization/MOM_shared_initialization.F90 index acffc9c927..2981bb9e94 100644 --- a/src/initialization/MOM_shared_initialization.F90 +++ b/src/initialization/MOM_shared_initialization.F90 @@ -11,7 +11,8 @@ module MOM_shared_initialization use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint use MOM_file_parser, only : get_param, log_param, param_file_type, log_version -use MOM_io, only : close_file, create_file, file_type, fieldtype, file_exists, field_size +use MOM_io, only : create_MOM_file, file_exists, field_size +use MOM_io, only : MOM_infra_file, MOM_field use MOM_io, only : MOM_read_data, MOM_read_vector, read_variable, stdout use MOM_io, only : open_file_to_read, close_file_to_read, SINGLE_FILE, MULTIPLE use MOM_io, only : slasher, vardesc, MOM_write_field, var_desc @@ -1346,9 +1347,9 @@ subroutine write_ocean_geometry_file(G, param_file, directory, US, geom_file) character(len=40) :: mdl = "write_ocean_geometry_file" type(vardesc), dimension(:), allocatable :: & vars ! Types with metadata about the variables and their staggering - type(fieldtype), dimension(:), allocatable :: & + type(MOM_field), dimension(:), allocatable :: & fields ! Opaque types used by MOM_io to store variable metadata information - type(file_type) :: IO_handle ! The I/O handle of the fileset + type(MOM_infra_file) :: IO_handle ! The I/O handle of the fileset integer :: nFlds ! The number of variables in this file integer :: file_threading logical :: multiple_files @@ -1412,7 +1413,8 @@ subroutine write_ocean_geometry_file(G, param_file, directory, US, geom_file) file_threading = SINGLE_FILE if (multiple_files) file_threading = MULTIPLE - call create_file(IO_handle, trim(filepath), vars, nFlds, fields, file_threading, dG=G) + call create_MOM_file(IO_handle, trim(filepath), vars, nFlds, fields, & + file_threading, dG=G) call MOM_write_field(IO_handle, fields(1), G%Domain, G%geoLatBu) call MOM_write_field(IO_handle, fields(2), G%Domain, G%geoLonBu) @@ -1445,7 +1447,7 @@ subroutine write_ocean_geometry_file(G, param_file, directory, US, geom_file) call MOM_write_field(IO_handle, fields(23), G%Domain, G%Dopen_v, scale=US%Z_to_m) endif - call close_file(IO_handle) + call IO_handle%close() deallocate(vars, fields)