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..45eb9ed065 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 + + +!> Testing a new approach +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..4604ae9e0b --- /dev/null +++ b/src/framework/MOM_io_file.F90 @@ -0,0 +1,1629 @@ +!> 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 + + +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 + + +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 + + +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 + + +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 + + +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 + + +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 + + +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 + + +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 + + +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 + + file_is_open_infra = fms2_file_is_open(handle%handle_infra) +end function file_is_open_infra + + +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 + type(MOM_field), intent(inout) :: fields(:) + + 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 + + +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 + + +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 + + +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 + + ! TODO: axis_nc is empty on non-root PEs! + 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 + + +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 + + ! TODO: field_nc is empty on non-root PEs! + 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 + + +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 + + +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 + + +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 + + +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 + + +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 + + +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 + + +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 + + file_is_open_nc = handle%is_open +end function file_is_open_nc + + +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 + + ! TODO: I could conditionally call and broadcast()... + 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 + type(MOM_field), intent(inout) :: fields(:) + + type(netcdf_axis), allocatable :: axes_nc(:) + type(netcdf_field), allocatable :: fields_nc(:) + integer :: i + + ! TODO: Support other ranks? + 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 + + +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 bde8e3e219..a592d92f0a 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 73be3f5843..39e0a2296a 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)