config_src/infra/FMS2/MOM_io_infra.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

1		! This file is part of MOM6, the Modular Ocean Model version 6.
2		! See the LICENSE file for licensing information.
3		! SPDX-License-Identifier: Apache-2.0
4
5		!> This module contains a thin inteface to mpp and fms I/O code
6		module MOM_io_infra
7
8		use MOM_domain_infra, only : MOM_domain_type, rescale_comp_data, AGRID, BGRID_NE, CGRID_NE
9		use MOM_domain_infra, only : domain2d, domain1d, CENTER, CORNER, NORTH_FACE, EAST_FACE
10		use MOM_error_infra, only : MOM_error=>MOM_err, NOTE, FATAL, WARNING, is_root_PE
11		use fms2_io_mod, only : fms2_open_file => open_file, check_if_open, fms2_close_file => close_file
12		use fms2_io_mod, only : fms2_flush_file => flush_file
13		use fms2_io_mod, only : FmsNetcdfDomainFile_t, FmsNetcdfFile_t, fms2_read_data => read_data
14		use fms2_io_mod, only : get_unlimited_dimension_name, get_num_dimensions, get_num_variables
15		use fms2_io_mod, only : get_variable_names, variable_exists, get_variable_size, get_variable_units
16		use fms2_io_mod, only : register_field, write_data, register_variable_attribute, register_global_attribute
17		use fms2_io_mod, only : variable_att_exists, get_variable_attribute, get_variable_num_dimensions
18		use fms2_io_mod, only : get_variable_dimension_names, is_dimension_registered, get_dimension_size
19		use fms2_io_mod, only : is_dimension_unlimited, register_axis, unlimited
20		use fms2_io_mod, only : get_dimension_names
21		use fms2_io_mod, only : get_global_io_domain_indices
22		use fms_io_utils_mod, only : fms2_file_exist => file_exists
23		use fms_io_utils_mod, only : get_filename_appendix
24
25		use fms_mod, only : write_version_number, check_nml_error
26		use mpp_domains_mod, only : mpp_get_compute_domain, mpp_get_global_domain
27		use mpp_mod, only : stdout_if_root=>stdout
28		use mpp_mod, only : mpp_pe, mpp_root_pe, mpp_npes
29		use mpp_mod, only : mpp_get_current_pelist_name
30		use iso_fortran_env, only : int64
31
32		implicit none ; private
33
34		! Duplication of FMS1 parameter values
35		! NOTE: Only kept to emulate FMS1 behavior, and may be removed in the future.
36		integer, parameter :: WRITEONLY_FILE = 100
37		integer, parameter :: READONLY_FILE = 101
38		integer, parameter :: APPEND_FILE = 102
39		integer, parameter :: OVERWRITE_FILE = 103
40		integer, parameter :: ASCII_FILE = 200
41		integer, parameter :: NETCDF_FILE = 203
42		integer, parameter :: SINGLE_FILE = 400
43		integer, parameter :: MULTIPLE = 401
44
45		! These interfaces are actually implemented or have explicit interfaces in this file.
46		public :: open_file, open_ASCII_file, file_is_open, close_file, flush_file, file_exists
47		public :: get_file_info, get_file_fields, get_file_times, get_filename_suffix
48		public :: read_field, read_vector, write_metadata, write_field
49		public :: field_exists, get_field_atts, get_field_size, read_field_chksum
50		public :: get_axis_size, get_axis_data, set_axis_data
51		public :: io_infra_init, io_infra_end, MOM_namelist_file, check_namelist_error, write_version
52		public :: stdout_if_root
53		! These types act as containers for information about files, fields and axes, respectively,
54		! and may also wrap opaque types from the underlying infrastructure.
55		public :: file_type, fieldtype, axistype
56		! These are encoding constant parmeters.
57		public :: ASCII_FILE, NETCDF_FILE, SINGLE_FILE, MULTIPLE
58		public :: APPEND_FILE, READONLY_FILE, OVERWRITE_FILE, WRITEONLY_FILE
59		public :: CENTER, CORNER, NORTH_FACE, EAST_FACE
60
61		!> Indicate whether a file exists, perhaps with domain decomposition
62		interface file_exists
63		module procedure FMS_file_exists
64		module procedure MOM_file_exists
65		end interface
66
67		!> Read a data field from a file
68		interface read_field
69		module procedure read_field_4d
70		module procedure read_field_3d, read_field_3d_region
71		module procedure read_field_2d, read_field_2d_region
72		module procedure read_field_1d, read_field_1d_int
73		module procedure read_field_0d, read_field_0d_int
74		end interface
75
76		!> Write a registered field to an output file
77		interface write_field
78		module procedure write_field_4d
79		module procedure write_field_3d
80		module procedure write_field_2d
81		module procedure write_field_1d
82		module procedure write_field_0d
83		module procedure MOM_write_axis
84		end interface write_field
85
86		!> Read a pair of data fields representing the two components of a vector from a file
87		interface read_vector
88		module procedure read_vector_3d
89		module procedure read_vector_2d
90		end interface read_vector
91
92		!> Write metadata about a variable or axis to a file and store it for later reuse
93		interface write_metadata
94		module procedure write_metadata_axis, write_metadata_field, write_metadata_global
95		end interface write_metadata
96
97		!> Close a file (or fileset). If the file handle does not point to an open file,
98		!! close_file simply returns without doing anything.
99		interface close_file
100		module procedure close_file_type, close_file_unit
101		end interface close_file
102
103		!> Type for holding a handle to an open file and related information
104		type :: file_type ; private
105		integer :: unit = -1 !< The framework identfier or netCDF unit number of an output file
106		type(FmsNetcdfDomainFile_t), pointer :: fileobj => NULL() !< A domain-decomposed
107		!! file object that is open for writing
108		character(len=:), allocatable :: filename !< The path to this file, if it is open
109		logical :: open_to_read = .false. !< If true, this file or fileset can be read
110		logical :: open_to_write = .false. !< If true, this file or fileset can be written to
111		integer :: num_times !< The number of time levels in this file
112		real :: file_time !< The time of the latest entry in the file.
113		end type file_type
114
115		!> This type is a container for information about a variable in a file.
116		type :: fieldtype ; private
117		character(len=256) :: name !< The name of this field in the files.
118		character(len=:), allocatable :: longname !< The long name for this field
119		character(len=:), allocatable :: units !< The units for this field
120		integer(kind=int64) :: chksum_read !< A checksum that has been read from a file
121		logical :: valid_chksum !< If true, this field has a valid checksum value.
122		end type fieldtype
123
124		!> This type is a container for information about an axis in a file.
125		type :: axistype ; private
126		character(len=256) :: name !< The name of this axis in the files.
127		real, allocatable, dimension(:) :: ax_data !< The values of the data on the axis.
128		logical :: domain_decomposed = .false. !< True if axis is domain-decomposed
129		end type axistype
130
131		contains
132
133		!> Reads the checksum value for a field that was recorded in a file, along with a flag indicating
134		!! whether the file contained a valid checksum for this field.
135	0	subroutine read_field_chksum(field, chksum, valid_chksum)
136		type(fieldtype), intent(in) :: field !< The field whose checksum attribute is to be read.
137		integer(kind=int64), intent(out) :: chksum !< The checksum for the field.
138		logical, intent(out) :: valid_chksum !< If true, chksum has been successfully read.
139
140	0	chksum = -1
141	0	valid_chksum = field%valid_chksum
142	0	if (valid_chksum) chksum = field%chksum_read
143
144	0	end subroutine read_field_chksum
145
146		!> Returns true if the named file or its domain-decomposed variant exists.
147	6	logical function MOM_file_exists(filename, MOM_Domain)
148		character(len=*), intent(in) :: filename !< The name of the file being inquired about
149		type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition
150
151	6	type(FmsNetcdfDomainFile_t) :: fileobj
152
153	12	MOM_file_exists = fms2_open_file(fileobj, filename, "read", MOM_Domain%mpp_domain)
154	6	if (MOM_file_exists) call fms2_close_file(fileobj)
155	12	end function MOM_file_exists
156
157		!> Returns true if the named file or its domain-decomposed variant exists.
158	6	logical function FMS_file_exists(filename)
159		character(len=*), intent(in) :: filename !< The name of the file being inquired about
160		! This function uses the fms_io function file_exist to determine whether
161		! a named file (or its decomposed variant) exists.
162
163	12	FMS_file_exists = fms2_file_exist(filename)
164	6	end function FMS_file_exists
165
166		!> indicates whether an I/O handle is attached to an open file
167	0	logical function file_is_open(IO_handle)
168		type(file_type), intent(in) :: IO_handle !< Handle to a file to inquire about
169
170	0	file_is_open = ((IO_handle%unit >= 0) .or. associated(IO_handle%fileobj))
171	0	end function file_is_open
172
173		!> closes a file (or fileset). If the file handle does not point to an open file,
174		!! close_file_type simply returns without doing anything.
175	2	subroutine close_file_type(IO_handle)
176		type(file_type), intent(inout) :: IO_handle !< The I/O handle for the file to be closed
177
178	2	if (associated(IO_handle%fileobj)) then
179	2	call fms2_close_file(IO_handle%fileobj)
180	2	deallocate(IO_handle%fileobj)
181		endif
182	2	if (allocated(IO_handle%filename)) deallocate(IO_handle%filename)
183	2	IO_handle%open_to_read = .false. ; IO_handle%open_to_write = .false.
184	2	IO_handle%num_times = 0 ; IO_handle%file_time = 0.0
185	2	end subroutine close_file_type
186
187		! TODO: close_file_unit is only used for ASCII files, which are opened outside
188		! of the framework, so this could probably be removed, and those calls could
189		! just be replaced with close(unit).
190
191		!> closes a file. If the unit does not point to an open file,
192		!! close_file_unit simply returns without doing anything.
193	6	subroutine close_file_unit(iounit)
194		integer, intent(inout) :: iounit !< The I/O unit for the file to be closed
195
196		logical :: unit_is_open
197
198	6	inquire(iounit, opened=unit_is_open)
199	6	if (unit_is_open) close(iounit)
200	6	end subroutine close_file_unit
201
202		!> Ensure that the output stream associated with a file handle is fully sent to disk.
203	0	subroutine flush_file(IO_handle)
204		type(file_type), intent(in) :: IO_handle !< The I/O handle for the file to flush
205
206	0	if (associated(IO_handle%fileobj)) then
207	0	call fms2_flush_file(IO_handle%fileobj)
208		endif
209	0	end subroutine flush_file
210
211		!> Initialize the underlying I/O infrastructure
212	1	subroutine io_infra_init(maxunits)
213		integer, optional, intent(in) :: maxunits !< An optional maximum number of file
214		!! unit numbers that can be used.
215
216		! FMS2 requires no explicit initialization, so this is a null function.
217	1	end subroutine io_infra_init
218
219		!> Gracefully close out and terminate the underlying I/O infrastructure
220	1	subroutine io_infra_end()
221		! FMS2 requires no explicit finalization, so this is a null function.
222	1	end subroutine io_infra_end
223
224		!> Open a single namelist file that is potentially readable by all PEs.
225	2	function MOM_namelist_file(filepath) result(iounit)
226		character(len=*), optional, intent(in) :: filepath
227		!< The file to open, by default "input.nml".
228		integer :: iounit
229		!< The opened unit number of the namelist file
230
231	2	character(len=:), allocatable :: nmlpath
232		! Namelist path
233	2	character(len=:), allocatable :: nmlpath_pe
234		! Hypothetical namelist path exclusive to the current PE list
235
236	2	if (present(filepath)) then
237	2	nmlpath = trim(filepath)
238		else
239		! FMS1 first checks for a namelist unique to the PE list, `input_{}.nml`.
240		! If not found, it defaults to `input.nml`.
241	0	nmlpath_pe = 'input_' // trim(mpp_get_current_pelist_name()) // '.nml'
242	0	if (file_exists(nmlpath_pe)) then
243	0	nmlpath = nmlpath_pe
244		else
245	0	nmlpath = 'input.nml'
246		endif
247		endif
248	2	call open_ASCII_file(iounit, nmlpath, action=READONLY_FILE)
249	4	end function MOM_namelist_file
250
251		!> Checks the iostat argument that is returned after reading a namelist variable and writes a
252		!! message if there is an error.
253	3	subroutine check_namelist_error(IOstat, nml_name)
254		integer, intent(in) :: IOstat !< An I/O status field from a namelist read call
255		character(len=*), intent(in) :: nml_name !< The name of the namelist
256		integer :: ierr
257	3	ierr = check_nml_error(IOstat, nml_name)
258	3	end subroutine check_namelist_error
259
260		!> Write a file version number to the log file or other output file
261	0	subroutine write_version(version, tag, unit)
262		character(len=*), intent(in) :: version !< A string that contains the routine name and version
263		character(len=*), optional, intent(in) :: tag !< A tag name to add to the message
264		integer, optional, intent(in) :: unit !< An alternate unit number for output
265
266	0	call write_version_number(version, tag, unit)
267	0	end subroutine write_version
268
269		!> open_file opens a file for parallel or single-file I/O.
270	2	subroutine open_file(IO_handle, filename, action, MOM_domain, threading, fileset)
271		type(file_type), intent(inout) :: IO_handle !< The handle for the opened file
272		character(len=*), intent(in) :: filename !< The path name of the file being opened
273		integer, optional, intent(in) :: action !< A flag indicating whether the file can be read
274		!! or written to and how to handle existing files.
275		!! The default is WRITE_ONLY.
276		type(MOM_domain_type), &
277		optional, intent(in) :: MOM_Domain !< A MOM_Domain that describes the decomposition
278		integer, optional, intent(in) :: threading !< A flag indicating whether one (SINGLE_FILE)
279		!! or multiple PEs (MULTIPLE) participate in I/O.
280		!! With the default, the root PE does I/O.
281		integer, optional, intent(in) :: fileset !< A flag indicating whether multiple PEs doing I/O due
282		!! to threading=MULTIPLE write to the same file (SINGLE_FILE)
283		!! or to one file per PE (MULTIPLE, the default).
284
285		! Local variables
286	2	type(FmsNetcdfDomainFile_t) :: fileobj_read ! A handle to a domain-decomposed file for obtaining information
287		! about the exiting time axis entries in append mode.
288		logical :: success ! If true, the file was opened successfully
289		integer :: file_mode ! An integer that encodes whether the file is to be opened for
290		! reading, writing or appending
291		character(len=40) :: mode ! A character string that encodes whether the file is to be opened for
292		! reading, writing or appending
293	2	character(len=:), allocatable :: filename_tmp ! A copy of filename with .nc appended if necessary.
294		character(len=256) :: dim_unlim_name ! name of the unlimited dimension in the file
295		integer :: index_nc
296
297	2	if (IO_handle%open_to_write) then
298		call MOM_error(WARNING, "open_file called for file "//trim(filename)//&
299	0	" with an IO_handle that is already open to to write.")
300	0	return
301		endif
302	2	if (IO_handle%open_to_read) then
303		call MOM_error(FATAL, "open_file called for file "//trim(filename)//&
304	0	" with an IO_handle that is already open to to read.")
305		endif
306
307	2	file_mode = WRITEONLY_FILE ; if (present(action)) file_mode = action
308
309		! Domains are currently required to use FMS I/O.
310		! NOTE: We restrict FMS2 IO usage to domain-based files due to issues with
311		! string-based attributes in certain compilers.
312		! But we may relax this requirement in the future.
313	2	if (.not. present(MOM_Domain)) &
314	0	call MOM_error(FATAL, 'open_file: FMS I/O requires a domain input.')
315
316	2	if (.not.associated(IO_handle%fileobj)) allocate (IO_handle%fileobj)
317
318		! The FMS1 interface automatically appends .nc if necessary, but FMS2 interface does not.
319	2	index_nc = index(trim(filename), ".nc")
320	2	if (index_nc > 0) then
321	2	filename_tmp = trim(filename)
322		else
323	0	filename_tmp = trim(filename)//".nc"
324	0	if (is_root_PE()) call MOM_error(WARNING, "Open_file is appending .nc to the filename "//trim(filename))
325		endif
326
327	2	if (file_mode == WRITEONLY_FILE) then ; mode = "write"
328	2	elseif (file_mode == APPEND_FILE) then ; mode = "append"
329	2	elseif (file_mode == OVERWRITE_FILE) then ; mode = "overwrite"
330	0	elseif (file_mode == READONLY_FILE) then ; mode = "read"
331		else
332	0	call MOM_error(FATAL, "open_file called with unrecognized action.")
333		endif
334
335	2	IO_handle%num_times = 0
336	2	IO_handle%file_time = 0.0
337	2	if ((file_mode == APPEND_FILE) .and. file_exists(filename_tmp, MOM_Domain)) then
338		! Determine the latest file time and number of records so far.
339	0	success = fms2_open_file(fileObj_read, trim(filename_tmp), "read", MOM_domain%mpp_domain)
340	0	dim_unlim_name = find_unlimited_dimension_name(fileObj_read)
341	0	if (len_trim(dim_unlim_name) > 0) &
342	0	call get_dimension_size(fileObj_read, trim(dim_unlim_name), IO_handle%num_times)
343	0	if (IO_handle%num_times > 0) &
344		call fms2_read_data(fileObj_read, trim(dim_unlim_name), IO_handle%file_time, &
345	0	unlim_dim_level=IO_handle%num_times)
346	0	call fms2_close_file(fileObj_read)
347		endif
348
349	2	success = fms2_open_file(IO_handle%fileobj, trim(filename_tmp), trim(mode), MOM_domain%mpp_domain)
350	2	if (.not.success) call MOM_error(FATAL, "Unable to open file "//trim(filename_tmp))
351	2	IO_handle%filename = trim(filename)
352
353	2	if (file_mode == READONLY_FILE) then
354	0	IO_handle%open_to_read = .true. ; IO_handle%open_to_write = .false.
355		else
356	2	IO_handle%open_to_read = .false. ; IO_handle%open_to_write = .true.
357		endif
358
359	4	end subroutine open_file
360
361		!> open_file opens an ascii file for parallel or single-file I/O using Fortran read and write calls.
362	8	subroutine open_ASCII_file(unit, file, action, threading, fileset)
363		integer, intent(out) :: unit !< The I/O unit for the opened file
364		character(len=*), intent(in) :: file !< The name of the file being opened
365		integer, optional, intent(in) :: action !< A flag indicating whether the file can be read
366		!! or written to and how to handle existing files.
367		integer, optional, intent(in) :: threading !< A flag indicating whether one (SINGLE_FILE)
368		!! or multiple PEs (MULTIPLE) participate in I/O.
369		!! With the default, the root PE does I/O.
370		integer, optional, intent(in) :: fileset !< A flag indicating whether multiple PEs doing I/O due
371		!! to threading=MULTIPLE write to the same file (SINGLE_FILE)
372		!! or to one file per PE (MULTIPLE, the default).
373
374		integer :: action_flag
375		integer :: threading_flag
376		integer :: fileset_flag
377		logical :: exists
378		logical :: is_open
379		character(len=6) :: action_arg, position_arg
380	8	character(len=:), allocatable :: filename
381
382		! NOTE: This function is written to emulate the original behavior of mpp_open
383		! from the FMS1 library, on which the MOM API is still based. Much of this
384		! can be removed if we choose to drop this compatibility, but for now we
385		! try to retain as much as possible.
386
387		! NOTE: Default FMS1 I/O settings are summarized below.
388		!
389		! access: Fortran and mpp_open default to SEQUENTIAL.
390		! form: The Fortran and mpp_open default (for MPP_ASCII) is FORMATTED.
391		! recl: mpp_open uses Fortran defaults when unset, so can be ignored.
392		! ios: FMS1 allowed this to be caught, but we do not support it.
393		! action/position: In mpp_open, these are inferred from `action`.
394		!
395		! MOM flag FMS1 flag action position
396		! -------- -------- ------ --------
397		! READONLY_FILE MPP_RDONLY READ REWIND
398		! WRITEONLY_FILE MPP_WRONLY WRITE REWIND
399		! OVERWRITE_FILE MPP_OVERWR WRITE REWIND
400		! APPEND_FILE MPP_APPEND WRITE APPEND
401		!
402		! From this, we can omit `access`, `form`, and `recl`, and can construct
403		! `action` and `position` from the input arguments.
404
405		! I/O configuration
406
407	8	action_flag = WRITEONLY_FILE
408	6	if (present(action)) action_flag = action
409
410	8	action_arg = 'write'
411	8	if (action_flag == READONLY_FILE) action_arg = 'read'
412
413	8	position_arg = 'rewind'
414	8	if (action_flag == APPEND_FILE) position_arg = 'append'
415
416		! Threading configuration
417
418	8	threading_flag = SINGLE_FILE
419	8	if (present(threading)) threading_flag = threading
420
421	8	fileset_flag = MULTIPLE
422	8	if (present(fileset)) fileset_flag = fileset
423
424		! Force fileset to be consistent with threading (as in FMS1)
425	8	if (threading_flag == SINGLE_FILE) fileset_flag = SINGLE_FILE
426
427		! Construct the distributed filename, if needed
428	8	filename = file
429	8	if (fileset_flag == MULTIPLE) then
430	0	if (mpp_npes() > 10000) then
431	0	write(filename, '(a,".",i6.6)') trim(filename), mpp_pe() - mpp_root_pe()
432		else
433	0	write(filename, '(a,".",i4.4)') trim(filename), mpp_pe() - mpp_root_pe()
434		endif
435		endif
436
437	8	inquire(file=filename, exist=exists)
438	8	if (exists .and. action_flag == WRITEONLY_FILE) &
439		call MOM_error(WARNING, 'open_ASCII_file: File ' // trim(filename) // &
440	0	' opened WRITEONLY already exists!')
441
442		open(newunit=unit, file=filename, action=trim(action_arg), &
443	8	position=trim(position_arg))
444
445		! This checks if open() failed but did not raise a runtime error.
446	8	inquire(unit, opened=is_open)
447	8	if (.not. is_open) &
448		call MOM_error(FATAL, &
449	0	'open_ASCII_file: File "' // trim(filename) // '" failed to open.')
450
451		! NOTE: There are two possible mpp_write_meta functions in FMS1:
452		! - call mpp_write_meta( unit, 'filename', cval=mpp_file(unit)%name)
453		! - call mpp_write_meta( unit, 'NumFilesInSet', ival=nfiles)
454		! I'm not convinced we actually want these, but note them here in case.
455	16	end subroutine open_ASCII_file
456
457
458		!> Provide a string to append to filenames, to differentiate ensemble members, for example.
459	7	subroutine get_filename_suffix(suffix)
460		character(len=*), intent(out) :: suffix !< A string to append to filenames
461
462	7	call get_filename_appendix(suffix)
463	7	end subroutine get_filename_suffix
464
465
466		!> Get information about the number of dimensions, variables and time levels
467		!! in the file associated with an open file unit
468	0	subroutine get_file_info(IO_handle, ndim, nvar, ntime)
469		type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for I/O
470		integer, optional, intent(out) :: ndim !< The number of dimensions in the file
471		integer, optional, intent(out) :: nvar !< The number of variables in the file
472		integer, optional, intent(out) :: ntime !< The number of time levels in the file
473
474		! Local variables
475		character(len=256) :: dim_unlim_name ! name of the unlimited dimension in the file
476		integer :: ndims, nvars, natts, ntimes
477
478	0	if (present(ndim)) ndim = get_num_dimensions(IO_handle%fileobj)
479	0	if (present(nvar)) nvar = get_num_variables(IO_handle%fileobj)
480	0	if (present(ntime)) then
481	0	ntime = 0
482	0	dim_unlim_name = find_unlimited_dimension_name(IO_handle%fileobj)
483	0	if (len_trim(dim_unlim_name) > 0) &
484	0	call get_dimension_size(IO_handle%fileobj, trim(dim_unlim_name), ntime)
485		endif
486	0	end subroutine get_file_info
487
488
489		!> Get the times of records from a file
490	0	subroutine get_file_times(IO_handle, time_values, ntime)
491		type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for I/O
492		real, allocatable, dimension(:), intent(inout) :: time_values !< The real times for the records in file.
493		integer, optional, intent(out) :: ntime !< The number of time levels in the file
494
495		character(len=256) :: dim_unlim_name ! name of the unlimited dimension in the file
496		integer :: ntimes ! The number of time levels in the file
497
498		!### Modify this routine to optionally convert to time_type, using information about the dimensions?
499
500	0	if (allocated(time_values)) deallocate(time_values)
501	0	call get_file_info(IO_handle, ntime=ntimes)
502	0	if (present(ntime)) ntime = ntimes
503	0	if (ntimes > 0) then
504	0	allocate(time_values(ntimes))
505	0	dim_unlim_name = find_unlimited_dimension_name(IO_handle%fileobj)
506	0	if (len_trim(dim_unlim_name) > 0) &
507	0	call fms2_read_data(IO_handle%fileobj, trim(dim_unlim_name), time_values)
508		endif
509	0	end subroutine get_file_times
510
511		!> Set up the field information (e.g., names and metadata) for all of the variables in a file. The
512		!! argument fields must be allocated with a size that matches the number of variables in a file.
513	0	subroutine get_file_fields(IO_handle, fields)
514		type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for I/O
515		type(fieldtype), dimension(:), intent(inout) :: fields !< Field-type descriptions of all of
516		!! the variables in a file.
517	0	character(len=256), dimension(size(fields)) :: var_names ! The names of all variables
518		character(len=256) :: units ! The units of a variable as recorded in the file
519		character(len=2048) :: longname ! The long-name of a variable as recorded in the file
520		character(len=64) :: checksum_char ! The hexadecimal checksum read from the file
521		integer(kind=int64), dimension(3) :: checksum_file ! The checksums for a variable in the file
522		integer :: nvar ! The number of variables in the file
523		integer :: i
524
525	0	nvar = size(fields)
526		! Local variables
527	0	call get_variable_names(IO_handle%fileobj, var_names)
528	0	do i=1,nvar
529	0	fields(i)%name = trim(var_names(i))
530	0	longname = ""
531	0	if (variable_att_exists(IO_handle%fileobj, var_names(i), "long_name")) &
532	0	call get_variable_attribute(IO_handle%fileobj, var_names(i), "long_name", longname)
533	0	fields(i)%longname = trim(longname)
534	0	units = ""
535	0	if (variable_att_exists(IO_handle%fileobj, var_names(i), "units")) &
536	0	call get_variable_attribute(IO_handle%fileobj, var_names(i), "units", units)
537	0	fields(i)%units = trim(units)
538
539	0	fields(i)%valid_chksum = variable_att_exists(IO_handle%fileobj, var_names(i), "checksum")
540	0	if (fields(i)%valid_chksum) then
541	0	call get_variable_attribute(IO_handle%fileobj, var_names(i), 'checksum', checksum_char)
542		! If there are problems, there might need to be code added to handle commas.
543	0	read (checksum_char(1:16), '(Z16)') fields(i)%chksum_read
544		endif
545		enddo
546	0	end subroutine get_file_fields
547
548		!> Extract information from a field type, as stored or as found in a file
549	0	subroutine get_field_atts(field, name, units, longname, checksum)
550		type(fieldtype), intent(in) :: field !< The field to extract information from
551		character(len=*), optional, intent(out) :: name !< The variable name
552		character(len=*), optional, intent(out) :: units !< The units of the variable
553		character(len=*), optional, intent(out) :: longname !< The long name of the variable
554		integer(kind=int64), dimension(:), &
555		optional, intent(out) :: checksum !< The checksums of the variable in a file
556
557	0	if (present(name)) name = trim(field%name)
558	0	if (present(units)) units = trim(field%units)
559	0	if (present(longname)) longname = trim(field%longname)
560	0	if (present(checksum)) checksum = field%chksum_read
561
562	0	end subroutine get_field_atts
563
564		!> Field_exists returns true if the field indicated by field_name is present in the
565		!! file file_name. If file_name does not exist, it returns false.
566	0	function field_exists(filename, field_name, domain, no_domain, MOM_domain)
567		character(len=*), intent(in) :: filename !< The name of the file being inquired about
568		character(len=*), intent(in) :: field_name !< The name of the field being sought
569		type(domain2d), target, optional, intent(in) :: domain !< A domain2d type that describes the decomposition
570		logical, optional, intent(in) :: no_domain !< This file does not use domain decomposition
571		type(MOM_domain_type), optional, intent(in) :: MOM_Domain !< A MOM_Domain that describes the decomposition
572		logical :: field_exists !< True if filename exists and field_name is in filename
573
574		! Local variables
575	0	type(FmsNetcdfDomainFile_t) :: fileObj_dd ! A handle to a domain-decomposed file for obtaining information
576		! about the exiting time axis entries in append mode.
577	0	type(FmsNetcdfFile_t) :: fileObj_simple ! A handle to a non-domain-decomposed file for obtaining information
578		! about the exiting time axis entries in append mode.
579		logical :: success ! If true, the file was opened successfully
580		logical :: domainless ! If true, this file does not use a domain-decomposed file.
581
582	0	domainless = .not.(present(MOM_domain) .or. present(domain))
583	0	if (present(no_domain)) then
584	0	if (domainless .and. .not.no_domain) call MOM_error(FATAL, &
585		"field_exists: When no_domain is present and false, a domain must be supplied in query about "//&
586	0	trim(field_name)//" in file "//trim(filename))
587	0	domainless = no_domain
588		endif
589
590	0	field_exists = .false.
591	0	if (file_exists(filename)) then
592	0	if (domainless) then
593	0	success = fms2_open_file(fileObj_simple, trim(filename), "read")
594	0	if (success) then
595	0	field_exists = variable_exists(fileObj_simple, field_name)
596	0	call fms2_close_file(fileObj_simple)
597		endif
598		else
599	0	if (present(MOM_domain)) then
600	0	success = fms2_open_file(fileObj_dd, trim(filename), "read", MOM_domain%mpp_domain)
601		else
602	0	success = fms2_open_file(fileObj_dd, trim(filename), "read", domain)
603		endif
604	0	if (success) then
605	0	field_exists = variable_exists(fileobj_dd, field_name)
606	0	call fms2_close_file(fileObj_dd)
607		endif
608		endif
609		endif
610	0	end function field_exists
611
612		!> Given filename and fieldname, this subroutine returns the size of the field in the file
613	0	subroutine get_field_size(filename, fieldname, sizes, field_found, no_domain)
614		character(len=*), intent(in) :: filename !< The name of the file to read
615		character(len=*), intent(in) :: fieldname !< The name of the variable whose sizes are returned
616		integer, dimension(:), intent(inout) :: sizes !< The sizes of the variable in each dimension
617		logical, optional, intent(out) :: field_found !< This indicates whether the field was found in
618		!! the input file. Without this argument, there
619		!! is a fatal error if the field is not found.
620		logical, optional, intent(in) :: no_domain !< If present and true, do not check for file
621		!! names with an appended tile number
622		! Local variables
623	0	type(FmsNetcdfFile_t) :: fileobj_read ! A handle to a non-domain-decomposed file for obtaining information
624		! about the exiting time axis entries in append mode.
625		logical :: success ! If true, the file was opened successfully
626		logical :: field_exists ! True if filename exists and field_name is in filename
627		integer :: i, ndims
628	0	character(len=512), allocatable :: dimnames(:) ! Field dimension names
629	0	logical, allocatable :: is_x(:), is_y(:), is_t(:) ! True if index matches axis type
630		integer :: size_indices(4) ! Mapping of size index to FMS1 convention
631		integer :: idx, swap
632
633	0	field_exists = .false.
634	0	if (file_exists(filename)) then
635	0	success = fms2_open_file(fileObj_read, trim(filename), "read")
636	0	if (success) then
637	0	field_exists = variable_exists(fileobj_read, fieldname)
638	0	if (field_exists) then
639	0	ndims = get_variable_num_dimensions(fileobj_read, fieldname)
640	0	if (ndims > size(sizes)) call MOM_error(FATAL, &
641	0	"get_field_size called with too few sizes for "//trim(fieldname)//" in "//trim(filename))
642	0	call get_variable_size(fileobj_read, fieldname, sizes(1:ndims))
643
644	0	do i=ndims+1,size(sizes) ; sizes(i) = 0 ; enddo
645
646		! If sizes exceeds ndims, then we fallback to the FMS1 convention
647		! where sizes has at least 4 dimension, and try to position values.
648	0	if (size(sizes) > ndims) then
649		! Assume FMS1 positioning rules: (nx, ny, nz, nt, ...)
650	0	if (size(sizes) < 4) &
651		call MOM_error(FATAL, "If sizes(:) exceeds field dimensions, "&
652	0	&"then its length must be at least 4.")
653
654		! Fall back to the FMS1 default values of 1 (from mpp field%size)
655	0	sizes(ndims+1:) = 1
656
657		! Gather the field dimension names
658	0	allocate(dimnames(ndims))
659	0	dimnames(:) = ""
660		call get_variable_dimension_names(fileObj_read, trim(fieldname), &
661	0	dimnames)
662
663		! Test the dimensions against standard (x,y,t) names and attributes
664	0	allocate(is_x(ndims), is_y(ndims), is_t(ndims))
665	0	is_x(:) = .false.
666	0	is_y(:) = .false.
667	0	is_t(:) = .false.
668		call categorize_axes(fileObj_read, filename, ndims, dimnames, &
669	0	is_x, is_y, is_t)
670
671		! Currently no z-test is supported, so disable assignment with 0
672		size_indices = [ &
673		find_index(is_x), &
674		find_index(is_y), &
675		0, &
676		find_index(is_t) &
677	0	]
678
679	0	do i = 1, size(size_indices)
680	0	idx = size_indices(i)
681	0	if (idx > 0) then
682	0	swap = sizes(i)
683	0	sizes(i) = sizes(idx)
684	0	sizes(idx) = swap
685		endif
686		enddo
687
688	0	deallocate(is_x, is_y, is_t)
689	0	deallocate(dimnames)
690		endif
691		endif
692		endif
693		endif
694	0	if (present(field_found)) field_found = field_exists
695	0	end subroutine get_field_size
696
697
698		!> Return the index of the first True element of a logical array.
699		!!
700		!! If all elements are false, return zero.
701	0	function find_index(vec) result(loc)
702		! NOTE: This function acts as a replacement for findloc() F2008 intrinsic,
703		! which is not available on some compilers, or may not support logicals.
704		logical, intent(in) :: vec(:)
705		integer :: loc
706
707		integer :: i
708
709	0	loc = 0
710	0	do i = 1, size(vec)
711	0	if (vec(i)) then
712	0	loc = i
713	0	exit
714		endif
715		enddo
716	0	end function find_index
717
718
719		!> Get the axis size from an axistype
720	0	function get_axis_size(axis) result(axis_size)
721		type(axistype), intent(in) :: axis
722		!< Infra axis
723		integer :: axis_size
724		!< Axis size
725
726	0	axis_size = size(axis%ax_data)
727	0	end function get_axis_size
728
729
730		!> Extracts and returns the axis data stored in an axistype.
731	0	subroutine get_axis_data(axis, axis_name, axis_data)
732		type(axistype), intent(in) :: axis
733		!< Infra axis
734		character(len=256), intent(out) :: axis_name
735		!< Axis name
736		real, dimension(:), intent(out) :: axis_data
737		!< Axis points
738
739		integer :: i
740
741	0	if (allocated(axis%ax_data)) then
742	0	if (size(axis%ax_data) > size(axis_data)) &
743		call MOM_error(FATAL, "get_axis_data called with too small of an " &
744	0	// "output data array for " // trim(axis%name) // ".")
745	0	do i=1,size(axis%ax_data)
746	0	axis_data(i) = axis%ax_data(i)
747		enddo
748		endif
749
750	0	axis_name = axis%name
751	0	end subroutine get_axis_data
752
753
754		!> Return a new axistype based on axis specs
755	0	subroutine set_axis_data(axis, axis_name, axis_data)
756		type(axistype), intent(inout) :: axis
757		!< Target axis
758		character(len=256), intent(in) :: axis_name
759		!< Target axis name
760		real, intent(in) :: axis_data(:)
761		!< Target axis values
762
763	0	axis%name = axis_name
764
765	0	if (allocated(axis%ax_data)) deallocate(axis%ax_data)
766	0	allocate(axis%ax_data(size(axis_data)))
767
768	0	axis%ax_data(:) = axis_data(:)
769
770		! NOTE: We do not yet consider domain-decomposed axes.
771	0	axis%domain_decomposed = .false.
772	0	end subroutine set_axis_data
773
774
775		!> This routine uses the fms_io subroutine read_data to read a scalar named
776		!! "fieldname" from a single or domain-decomposed file "filename".
777	0	subroutine read_field_0d(filename, fieldname, data, timelevel, scale, MOM_Domain, &
778		global_file, file_may_be_4d)
779		character(len=*), intent(in) :: filename !< The name of the file to read
780		character(len=*), intent(in) :: fieldname !< The variable name of the data in the file
781		real, intent(inout) :: data !< The 1-dimensional array into which the data
782		integer, optional, intent(in) :: timelevel !< The time level in the file to read
783		real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied
784		!! by before it is returned.
785		type(MOM_domain_type), &
786		optional, intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition
787		logical, optional, intent(in) :: global_file !< If true, read from a single global file
788		logical, optional, intent(in) :: file_may_be_4d !< If true, this file may have 4-d arrays, but
789		!! with the FMS2 I/O interfaces this does not matter.
790
791		! Local variables
792	0	type(FmsNetcdfFile_t) :: fileObj ! A handle to a non-domain-decomposed file
793	0	type(FmsNetcdfDomainFile_t) :: fileobj_DD ! A handle to a domain-decomposed file object
794		character(len=96) :: var_to_read ! Name of variable to read from the netcdf file
795		logical :: has_time_dim ! True if the variable has an unlimited time axis.
796		logical :: success ! True if the file was successfully opened
797
798	0	if (present(MOM_Domain)) then
799		! Open the FMS2 file-set.
800	0	success = fms2_open_file(fileobj_DD, filename, "read", MOM_domain%mpp_domain)
801	0	if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename))
802
803		! Find the matching case-insensitive variable name in the file and prepare to read it.
804		call prepare_to_read_var(fileobj_DD, fieldname, "read_field_0d: ", filename, &
805	0	var_to_read, has_time_dim, timelevel)
806
807		! Read the data.
808	0	if (present(timelevel) .and. has_time_dim) then
809	0	call fms2_read_data(fileobj_DD, var_to_read, data, unlim_dim_level=timelevel)
810		else
811	0	call fms2_read_data(fileobj_DD, var_to_read, data)
812		endif
813
814		! Close the file-set.
815	0	if (check_if_open(fileobj_DD)) call fms2_close_file(fileobj_DD)
816		else
817		! Open the FMS2 file-set.
818	0	success = fms2_open_file(fileObj, trim(filename), "read")
819	0	if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename))
820
821		! Find the matching case-insensitive variable name in the file, and determine whether it
822		! has a time dimension.
823		call find_varname_in_file(fileObj, fieldname, "read_field_0d: ", filename, &
824	0	var_to_read, has_time_dim, timelevel)
825
826		! Read the data.
827	0	if (present(timelevel) .and. has_time_dim) then
828	0	call fms2_read_data(fileobj, var_to_read, data, unlim_dim_level=timelevel)
829		else
830	0	call fms2_read_data(fileobj, var_to_read, data)
831		endif
832
833		! Close the file-set.
834	0	if (check_if_open(fileobj)) call fms2_close_file(fileobj)
835		endif
836
837	0	if (present(scale)) then ; if (scale /= 1.0) then
838	0	data = scale*data
839		endif ; endif
840
841	0	end subroutine read_field_0d
842
843		!> This routine uses the fms_io subroutine read_data to read a 1-D data field named
844		!! "fieldname" from a single or domain-decomposed file "filename".
845	0	subroutine read_field_1d(filename, fieldname, data, timelevel, scale, MOM_Domain, &
846		global_file, file_may_be_4d)
847		character(len=*), intent(in) :: filename !< The name of the file to read
848		character(len=*), intent(in) :: fieldname !< The variable name of the data in the file
849		real, dimension(:), intent(inout) :: data !< The 1-dimensional array into which the data
850		integer, optional, intent(in) :: timelevel !< The time level in the file to read
851		real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied
852		!! by before they are returned.
853		type(MOM_domain_type), &
854		optional, intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition
855		logical, optional, intent(in) :: global_file !< If true, read from a single global file
856		logical, optional, intent(in) :: file_may_be_4d !< If true, this file may have 4-d arrays, but
857		!! with the FMS2 I/O interfaces this does not matter.
858
859		! Local variables
860	0	type(FmsNetcdfFile_t) :: fileObj ! A handle to a non-domain-decomposed file
861	0	type(FmsNetcdfDomainFile_t) :: fileobj_DD ! A handle to a domain-decomposed file object
862		character(len=96) :: var_to_read ! Name of variable to read from the netcdf file
863		logical :: has_time_dim ! True if the variable has an unlimited time axis.
864		logical :: success ! True if the file was successfully opened
865
866	0	if (present(MOM_Domain)) then
867		! Open the FMS2 file-set.
868	0	success = fms2_open_file(fileobj_DD, filename, "read", MOM_domain%mpp_domain)
869	0	if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename))
870
871		! Find the matching case-insensitive variable name in the file and prepare to read it.
872		call prepare_to_read_var(fileobj_DD, fieldname, "read_field_1d: ", filename, &
873	0	var_to_read, has_time_dim, timelevel)
874
875		! Read the data.
876	0	if (present(timelevel) .and. has_time_dim) then
877	0	call fms2_read_data(fileobj_DD, var_to_read, data, unlim_dim_level=timelevel)
878		else
879	0	call fms2_read_data(fileobj_DD, var_to_read, data)
880		endif
881
882		! Close the file-set.
883	0	if (check_if_open(fileobj_DD)) call fms2_close_file(fileobj_DD)
884		else
885		! Open the FMS2 file-set.
886	0	success = fms2_open_file(fileObj, trim(filename), "read")
887	0	if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename))
888
889		! Find the matching case-insensitive variable name in the file, and determine whether it
890		! has a time dimension.
891		call find_varname_in_file(fileObj, fieldname, "read_field_1d: ", filename, &
892	0	var_to_read, has_time_dim, timelevel)
893
894		! Read the data.
895	0	if (present(timelevel) .and. has_time_dim) then
896	0	call fms2_read_data(fileobj, var_to_read, data, unlim_dim_level=timelevel)
897		else
898	0	call fms2_read_data(fileobj, var_to_read, data)
899		endif
900
901		! Close the file-set.
902	0	if (check_if_open(fileobj)) call fms2_close_file(fileobj)
903		endif
904
905	0	if (present(scale)) then ; if (scale /= 1.0) then
906	0	data(:) = scale*data(:)
907		endif ; endif
908
909	0	end subroutine read_field_1d
910
911		!> This routine uses the fms_io subroutine read_data to read a distributed
912		!! 2-D data field named "fieldname" from file "filename". Valid values for
913		!! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE.
914	0	subroutine read_field_2d(filename, fieldname, data, MOM_Domain, &
915		timelevel, position, scale, global_file, file_may_be_4d)
916		character(len=*), intent(in) :: filename !< The name of the file to read
917		character(len=*), intent(in) :: fieldname !< The variable name of the data in the file
918		real, dimension(:,:), intent(inout) :: data !< The 2-dimensional array into which the data
919		!! should be read
920		type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition
921		integer, optional, intent(in) :: timelevel !< The time level in the file to read
922		integer, optional, intent(in) :: position !< A flag indicating where this data is located
923		real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied
924		!! by before it is returned.
925		logical, optional, intent(in) :: global_file !< If true, read from a single global file
926		logical, optional, intent(in) :: file_may_be_4d !< If true, this file may have 4-d arrays, but
927		!! with the FMS2 I/O interfaces this does not matter.
928
929		! Local variables
930	0	type(FmsNetcdfDomainFile_t) :: fileobj ! A handle to a domain-decomposed file object
931		character(len=96) :: var_to_read ! Name of variable to read from the netcdf file
932		logical :: has_time_dim ! True if the variable has an unlimited time axis.
933		logical :: success ! True if the file was successfully opened
934
935		! Open the FMS2 file-set.
936	0	success = fms2_open_file(fileobj, filename, "read", MOM_domain%mpp_domain)
937	0	if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename))
938
939		! Find the matching case-insensitive variable name in the file and prepare to read it.
940		call prepare_to_read_var(fileobj, fieldname, "read_field_2d: ", filename, &
941	0	var_to_read, has_time_dim, timelevel, position)
942
943		! Read the data.
944	0	if (present(timelevel) .and. has_time_dim) then
945	0	call fms2_read_data(fileobj, var_to_read, data, unlim_dim_level=timelevel)
946		else
947	0	call fms2_read_data(fileobj, var_to_read, data)
948		endif
949
950		! Close the file-set.
951	0	if (check_if_open(fileobj)) call fms2_close_file(fileobj)
952
953	0	if (present(scale)) then ; if (scale /= 1.0) then
954	0	call rescale_comp_data(MOM_Domain, data, scale)
955		endif ; endif
956
957	0	end subroutine read_field_2d
958
959		!> This routine uses the fms_io subroutine read_data to read a region from a distributed or
960		!! global 2-D data field named "fieldname" from file "filename".
961	0	subroutine read_field_2d_region(filename, fieldname, data, start, nread, MOM_domain, &
962		no_domain, scale)
963		character(len=*), intent(in) :: filename !< The name of the file to read
964		character(len=*), intent(in) :: fieldname !< The variable name of the data in the file
965		real, dimension(:,:), intent(inout) :: data !< The 2-dimensional array into which the data
966		!! should be read
967		integer, dimension(:), intent(in) :: start !< The starting index to read in each of 4
968		!! dimensions. For this 2-d read, the 3rd
969		!! and 4th values are always 1.
970		integer, dimension(:), intent(in) :: nread !< The number of points to read in each of 4
971		!! dimensions. For this 2-d read, the 3rd
972		!! and 4th values are always 1.
973		type(MOM_domain_type), &
974		optional, intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition
975		logical, optional, intent(in) :: no_domain !< If present and true, this variable does not
976		!! use domain decomposion.
977		real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied
978		!! by before it is returned.
979
980		! Local variables
981	0	type(FmsNetcdfFile_t) :: fileObj ! A handle to a non-domain-decomposed file
982	0	type(FmsNetcdfDomainFile_t) :: fileobj_DD ! A handle to a domain-decomposed file object
983		character(len=96) :: var_to_read ! Name of variable to read from the netcdf file
984		logical :: success ! True if the file was successfully opened
985
986	0	if (present(MOM_Domain)) then
987		! Open the FMS2 file-set.
988	0	success = fms2_open_file(fileobj_DD, filename, "read", MOM_domain%mpp_domain)
989	0	if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename))
990
991		! Find the matching case-insensitive variable name in the file and prepare to read it.
992		call prepare_to_read_var(fileobj_DD, fieldname, "read_field_2d_region: ", &
993	0	filename, var_to_read)
994
995		! Read the data.
996	0	call fms2_read_data(fileobj_DD, var_to_read, data, corner=start(1:2), edge_lengths=nread(1:2))
997
998		! Close the file-set.
999	0	if (check_if_open(fileobj_DD)) call fms2_close_file(fileobj_DD)
1000		else
1001		! Open the FMS2 file-set.
1002	0	success = fms2_open_file(fileObj, trim(filename), "read")
1003	0	if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename))
1004
1005		! Find the matching case-insensitive variable name in the file, and determine whether it
1006		! has a time dimension.
1007	0	call find_varname_in_file(fileObj, fieldname, "read_field_2d_region: ", filename, var_to_read)
1008
1009		! Read the data.
1010	0	call fms2_read_data(fileobj, var_to_read, data, corner=start(1:2), edge_lengths=nread(1:2))
1011
1012		! Close the file-set.
1013	0	if (check_if_open(fileobj)) call fms2_close_file(fileobj)
1014		endif
1015
1016	0	if (present(scale)) then ; if (scale /= 1.0) then
1017	0	if (present(MOM_Domain)) then
1018	0	call rescale_comp_data(MOM_Domain, data, scale)
1019		else
1020		! Dangerously rescale the whole array
1021	0	data(:,:) = scale*data(:,:)
1022		endif
1023		endif ; endif
1024
1025	0	end subroutine read_field_2d_region
1026
1027		!> This routine uses the fms_io subroutine read_data to read a distributed
1028		!! 3-D data field named "fieldname" from file "filename". Valid values for
1029		!! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE.
1030	0	subroutine read_field_3d(filename, fieldname, data, MOM_Domain, &
1031		timelevel, position, scale, global_file, file_may_be_4d)
1032		character(len=*), intent(in) :: filename !< The name of the file to read
1033		character(len=*), intent(in) :: fieldname !< The variable name of the data in the file
1034		real, dimension(:,:,:), intent(inout) :: data !< The 3-dimensional array into which the data
1035		!! should be read
1036		type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition
1037		integer, optional, intent(in) :: timelevel !< The time level in the file to read
1038		integer, optional, intent(in) :: position !< A flag indicating where this data is located
1039		real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied
1040		!! by before it is returned.
1041		logical, optional, intent(in) :: global_file !< If true, read from a single global file
1042		logical, optional, intent(in) :: file_may_be_4d !< If true, this file may have 4-d arrays, but
1043		!! with the FMS2 I/O interfaces this does not matter.
1044
1045		! Local variables
1046	0	type(FmsNetcdfDomainFile_t) :: fileobj ! A handle to a domain-decomposed file object
1047		character(len=96) :: var_to_read ! Name of variable to read from the netcdf file
1048		logical :: has_time_dim ! True if the variable has an unlimited time axis.
1049		logical :: success ! True if the file was successfully opened
1050
1051		! Open the FMS2 file-set.
1052	0	success = fms2_open_file(fileobj, filename, "read", MOM_domain%mpp_domain)
1053	0	if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename))
1054
1055		! Find the matching case-insensitive variable name in the file and prepare to read it.
1056		call prepare_to_read_var(fileobj, fieldname, "read_field_3d: ", filename, &
1057	0	var_to_read, has_time_dim, timelevel, position)
1058
1059		! Read the data.
1060	0	if (present(timelevel) .and. has_time_dim) then
1061	0	call fms2_read_data(fileobj, var_to_read, data, unlim_dim_level=timelevel)
1062		else
1063	0	call fms2_read_data(fileobj, var_to_read, data)
1064		endif
1065
1066		! Close the file-set.
1067	0	if (check_if_open(fileobj)) call fms2_close_file(fileobj)
1068
1069	0	if (present(scale)) then ; if (scale /= 1.0) then
1070	0	call rescale_comp_data(MOM_Domain, data, scale)
1071		endif ; endif
1072
1073	0	end subroutine read_field_3d
1074
1075		!> This routine uses the fms_io subroutine read_data to read a region from a distributed or
1076		!! global 3-D data field named "fieldname" from file "filename".
1077	0	subroutine read_field_3d_region(filename, fieldname, data, start, nread, MOM_domain, &
1078		no_domain, scale)
1079		character(len=*), intent(in) :: filename !< The name of the file to read
1080		character(len=*), intent(in) :: fieldname !< The variable name of the data in the file
1081		real, dimension(:,:,:), intent(inout) :: data !< The 3-dimensional array into which the data
1082		!! should be read
1083		integer, dimension(:), intent(in) :: start !< The starting index to read in each of 3
1084		!! dimensions. For this 3-d read, the
1085		!! 4th value is always 1.
1086		integer, dimension(:), intent(in) :: nread !< The number of points to read in each of 4
1087		!! dimensions. For this 3-d read, the
1088		!! 4th values are always 1.
1089		type(MOM_domain_type), &
1090		optional, intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition
1091		logical, optional, intent(in) :: no_domain !< If present and true, this variable does not
1092		!! use domain decomposion.
1093		real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied
1094		!! by before it is returned.
1095
1096		! Local variables
1097	0	type(FmsNetcdfFile_t) :: fileObj ! A handle to a non-domain-decomposed file
1098	0	type(FmsNetcdfDomainFile_t) :: fileobj_DD ! A handle to a domain-decomposed file object
1099		character(len=96) :: var_to_read ! Name of variable to read from the netcdf file
1100		logical :: success ! True if the file was successfully opened
1101
1102	0	if (present(MOM_Domain)) then
1103		! Open the FMS2 file-set.
1104	0	success = fms2_open_file(fileobj_DD, filename, "read", MOM_domain%mpp_domain)
1105	0	if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename))
1106
1107		! Find the matching case-insensitive variable name in the file and prepare to read it.
1108		call prepare_to_read_var(fileobj_DD, fieldname, "read_field_2d_region: ", &
1109	0	filename, var_to_read)
1110
1111		! Read the data.
1112	0	call fms2_read_data(fileobj_DD, var_to_read, data, corner=start(1:3), edge_lengths=nread(1:3))
1113
1114		! Close the file-set.
1115	0	if (check_if_open(fileobj_DD)) call fms2_close_file(fileobj_DD)
1116		else
1117		! Open the FMS2 file-set.
1118	0	success = fms2_open_file(fileObj, trim(filename), "read")
1119	0	if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename))
1120
1121		! Find the matching case-insensitive variable name in the file, and determine whether it
1122		! has a time dimension.
1123	0	call find_varname_in_file(fileObj, fieldname, "read_field_2d_region: ", filename, var_to_read)
1124
1125		! Read the data.
1126	0	call fms2_read_data(fileobj, var_to_read, data, corner=start(1:3), edge_lengths=nread(1:3))
1127
1128		! Close the file-set.
1129	0	if (check_if_open(fileobj)) call fms2_close_file(fileobj)
1130		endif
1131
1132	0	if (present(scale)) then ; if (scale /= 1.0) then
1133	0	if (present(MOM_Domain)) then
1134	0	call rescale_comp_data(MOM_Domain, data, scale)
1135		else
1136		! Dangerously rescale the whole array
1137	0	data(:,:,:) = scale*data(:,:,:)
1138		endif
1139		endif ; endif
1140
1141	0	end subroutine read_field_3d_region
1142
1143		!> This routine uses the fms_io subroutine read_data to read a distributed
1144		!! 4-D data field named "fieldname" from file "filename". Valid values for
1145		!! "position" include CORNER, CENTER, EAST_FACE and NORTH_FACE.
1146	0	subroutine read_field_4d(filename, fieldname, data, MOM_Domain, &
1147		timelevel, position, scale, global_file)
1148		character(len=*), intent(in) :: filename !< The name of the file to read
1149		character(len=*), intent(in) :: fieldname !< The variable name of the data in the file
1150		real, dimension(:,:,:,:), intent(inout) :: data !< The 4-dimensional array into which the data
1151		!! should be read
1152		type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition
1153		integer, optional, intent(in) :: timelevel !< The time level in the file to read
1154		integer, optional, intent(in) :: position !< A flag indicating where this data is located
1155		real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied
1156		!! by before it is returned.
1157		logical, optional, intent(in) :: global_file !< If true, read from a single global file
1158
1159
1160		! Local variables
1161	0	type(FmsNetcdfDomainFile_t) :: fileobj ! A handle to a domain-decomposed file object
1162		logical :: has_time_dim ! True if the variable has an unlimited time axis.
1163		character(len=96) :: var_to_read ! Name of variable to read from the netcdf file
1164		logical :: success ! True if the file was successfully opened
1165
1166		! Open the FMS2 file-set.
1167	0	success = fms2_open_file(fileobj, filename, "read", MOM_domain%mpp_domain)
1168	0	if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename))
1169
1170		! Find the matching case-insensitive variable name in the file and prepare to read it.
1171		call prepare_to_read_var(fileobj, fieldname, "read_field_4d: ", filename, &
1172	0	var_to_read, has_time_dim, timelevel, position)
1173
1174		! Read the data.
1175	0	if (present(timelevel) .and. has_time_dim) then
1176	0	call fms2_read_data(fileobj, var_to_read, data, unlim_dim_level=timelevel)
1177		else
1178	0	call fms2_read_data(fileobj, var_to_read, data)
1179		endif
1180
1181		! Close the file-set.
1182	0	if (check_if_open(fileobj)) call fms2_close_file(fileobj)
1183
1184	0	if (present(scale)) then ; if (scale /= 1.0) then
1185	0	call rescale_comp_data(MOM_Domain, data, scale)
1186		endif ; endif
1187
1188	0	end subroutine read_field_4d
1189
1190		!> This routine uses the fms_io subroutine read_data to read a scalar integer
1191		!! data field named "fieldname" from file "filename".
1192	0	subroutine read_field_0d_int(filename, fieldname, data, timelevel)
1193		character(len=*), intent(in) :: filename !< The name of the file to read
1194		character(len=*), intent(in) :: fieldname !< The variable name of the data in the file
1195		integer, intent(inout) :: data !< The 1-dimensional array into which the data
1196		integer, optional, intent(in) :: timelevel !< The time level in the file to read
1197
1198		! Local variables
1199	0	type(FmsNetcdfFile_t) :: fileObj ! A handle to a non-domain-decomposed file
1200		logical :: has_time_dim ! True if the variable has an unlimited time axis.
1201		character(len=96) :: var_to_read ! Name of variable to read from the netcdf file
1202		logical :: success ! If true, the file was opened successfully
1203
1204		! This routine might not be needed for MOM6.
1205
1206		! Open the FMS2 file-set.
1207	0	success = fms2_open_file(fileObj, trim(filename), "read")
1208	0	if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename))
1209
1210		! Find the matching case-insensitive variable name in the file, and determine whether it
1211		! has a time dimension.
1212		call find_varname_in_file(fileObj, fieldname, "read_field_0d_int: ", filename, &
1213	0	var_to_read, has_time_dim, timelevel)
1214
1215		! Read the data.
1216	0	if (present(timelevel) .and. has_time_dim) then
1217	0	call fms2_read_data(fileobj, var_to_read, data, unlim_dim_level=timelevel)
1218		else
1219	0	call fms2_read_data(fileobj, var_to_read, data)
1220		endif
1221
1222		! Close the file-set.
1223	0	if (check_if_open(fileobj)) call fms2_close_file(fileobj)
1224	0	end subroutine read_field_0d_int
1225
1226		!> This routine uses the fms_io subroutine read_data to read a 1-D integer
1227		!! data field named "fieldname" from file "filename".
1228	0	subroutine read_field_1d_int(filename, fieldname, data, timelevel)
1229		character(len=*), intent(in) :: filename !< The name of the file to read
1230		character(len=*), intent(in) :: fieldname !< The variable name of the data in the file
1231		integer, dimension(:), intent(inout) :: data !< The 1-dimensional array into which the data
1232		integer, optional, intent(in) :: timelevel !< The time level in the file to read
1233
1234		! Local variables
1235	0	type(FmsNetcdfFile_t) :: fileObj ! A handle to a non-domain-decomposed file for obtaining information
1236		! about the exiting time axis entries in append mode.
1237		logical :: has_time_dim ! True if the variable has an unlimited time axis.
1238		character(len=96) :: var_to_read ! Name of variable to read from the netcdf file
1239		logical :: success ! If true, the file was opened successfully
1240
1241		! This routine might not be needed for MOM6.
1242
1243		! Open the FMS2 file-set.
1244	0	success = fms2_open_file(fileObj, trim(filename), "read")
1245	0	if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename))
1246
1247		! Find the matching case-insensitive variable name in the file, and determine whether it
1248		! has a time dimension.
1249		call find_varname_in_file(fileObj, fieldname, "read_field_1d_int: ", filename, &
1250	0	var_to_read, has_time_dim, timelevel)
1251
1252		! Read the data.
1253	0	if (present(timelevel) .and. has_time_dim) then
1254	0	call fms2_read_data(fileobj, var_to_read, data, unlim_dim_level=timelevel)
1255		else
1256	0	call fms2_read_data(fileobj, var_to_read, data)
1257		endif
1258
1259		! Close the file-set.
1260	0	if (check_if_open(fileobj)) call fms2_close_file(fileobj)
1261	0	end subroutine read_field_1d_int
1262
1263
1264		!> This routine uses the fms_io subroutine read_data to read a pair of distributed
1265		!! 2-D data fields with names given by "[uv]_fieldname" from file "filename". Valid values for
1266		!! "stagger" include CGRID_NE, BGRID_NE, and AGRID.
1267	0	subroutine read_vector_2d(filename, u_fieldname, v_fieldname, u_data, v_data, MOM_Domain, &
1268		timelevel, stagger, scalar_pair, scale)
1269		character(len=*), intent(in) :: filename !< The name of the file to read
1270		character(len=*), intent(in) :: u_fieldname !< The variable name of the u data in the file
1271		character(len=*), intent(in) :: v_fieldname !< The variable name of the v data in the file
1272		real, dimension(:,:), intent(inout) :: u_data !< The 2 dimensional array into which the
1273		!! u-component of the data should be read
1274		real, dimension(:,:), intent(inout) :: v_data !< The 2 dimensional array into which the
1275		!! v-component of the data should be read
1276		type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition
1277		integer, optional, intent(in) :: timelevel !< The time level in the file to read
1278		integer, optional, intent(in) :: stagger !< A flag indicating where this vector is discretized
1279		logical, optional, intent(in) :: scalar_pair !< If true, a pair of scalars are to be read
1280		real, optional, intent(in) :: scale !< A scaling factor that the fields are multiplied
1281		!! by before they are returned.
1282		! Local variables
1283	0	type(FmsNetcdfDomainFile_t) :: fileobj ! A handle to a domain-decomposed file object
1284		logical :: has_time_dim ! True if the variables have an unlimited time axis.
1285		character(len=96) :: u_var, v_var ! Name of u and v variables to read from the netcdf file
1286		logical :: success ! True if the file was successfully opened
1287		integer :: u_pos, v_pos ! Flags indicating the positions of the u- and v- components.
1288
1289	0	u_pos = EAST_FACE ; v_pos = NORTH_FACE
1290	0	if (present(stagger)) then
1291	0	if (stagger == CGRID_NE) then ; u_pos = EAST_FACE ; v_pos = NORTH_FACE
1292	0	elseif (stagger == BGRID_NE) then ; u_pos = CORNER ; v_pos = CORNER
1293	0	elseif (stagger == AGRID) then ; u_pos = CENTER ; v_pos = CENTER ; endif
1294		endif
1295
1296		! Open the FMS2 file-set.
1297	0	success = fms2_open_file(fileobj, filename, "read", MOM_domain%mpp_domain)
1298	0	if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename))
1299
1300		! Find the matching case-insensitive u- and v-variable names in the file and prepare to read them.
1301		call prepare_to_read_var(fileobj, u_fieldname, "read_vector_2d: ", filename, &
1302	0	u_var, has_time_dim, timelevel, position=u_pos)
1303		call prepare_to_read_var(fileobj, v_fieldname, "read_vector_2d: ", filename, &
1304	0	v_var, has_time_dim, timelevel, position=v_pos)
1305
1306		! Read the u-data and v-data. There would already been an error message for one
1307		! of the variables if they are inconsistent in having an unlimited dimension.
1308	0	if (present(timelevel) .and. has_time_dim) then
1309	0	call fms2_read_data(fileobj, u_var, u_data, unlim_dim_level=timelevel)
1310	0	call fms2_read_data(fileobj, v_var, v_data, unlim_dim_level=timelevel)
1311		else
1312	0	call fms2_read_data(fileobj, u_var, u_data)
1313	0	call fms2_read_data(fileobj, v_var, v_data)
1314		endif
1315
1316		! Close the file-set.
1317	0	if (check_if_open(fileobj)) call fms2_close_file(fileobj)
1318
1319	0	if (present(scale)) then ; if (scale /= 1.0) then
1320	0	call rescale_comp_data(MOM_Domain, u_data, scale)
1321	0	call rescale_comp_data(MOM_Domain, v_data, scale)
1322		endif ; endif
1323
1324	0	end subroutine read_vector_2d
1325
1326		!> This routine uses the fms_io subroutine read_data to read a pair of distributed
1327		!! 3-D data fields with names given by "[uv]_fieldname" from file "filename". Valid values for
1328		!! "stagger" include CGRID_NE, BGRID_NE, and AGRID.
1329	0	subroutine read_vector_3d(filename, u_fieldname, v_fieldname, u_data, v_data, MOM_Domain, &
1330		timelevel, stagger, scalar_pair, scale)
1331		character(len=*), intent(in) :: filename !< The name of the file to read
1332		character(len=*), intent(in) :: u_fieldname !< The variable name of the u data in the file
1333		character(len=*), intent(in) :: v_fieldname !< The variable name of the v data in the file
1334		real, dimension(:,:,:), intent(inout) :: u_data !< The 3 dimensional array into which the
1335		!! u-component of the data should be read
1336		real, dimension(:,:,:), intent(inout) :: v_data !< The 3 dimensional array into which the
1337		!! v-component of the data should be read
1338		type(MOM_domain_type), intent(in) :: MOM_Domain !< The MOM_Domain that describes the decomposition
1339		integer, optional, intent(in) :: timelevel !< The time level in the file to read
1340		integer, optional, intent(in) :: stagger !< A flag indicating where this vector is discretized
1341		logical, optional, intent(in) :: scalar_pair !< If true, a pair of scalars are to be read.
1342		real, optional, intent(in) :: scale !< A scaling factor that the fields are multiplied
1343		!! by before they are returned.
1344
1345		! Local variables
1346	0	type(FmsNetcdfDomainFile_t) :: fileobj ! A handle to a domain-decomposed file object
1347		logical :: has_time_dim ! True if the variables have an unlimited time axis.
1348		character(len=96) :: u_var, v_var ! Name of u and v variables to read from the netcdf file
1349		logical :: success ! True if the file was successfully opened
1350		integer :: u_pos, v_pos ! Flags indicating the positions of the u- and v- components.
1351
1352	0	u_pos = EAST_FACE ; v_pos = NORTH_FACE
1353	0	if (present(stagger)) then
1354	0	if (stagger == CGRID_NE) then ; u_pos = EAST_FACE ; v_pos = NORTH_FACE
1355	0	elseif (stagger == BGRID_NE) then ; u_pos = CORNER ; v_pos = CORNER
1356	0	elseif (stagger == AGRID) then ; u_pos = CENTER ; v_pos = CENTER ; endif
1357		endif
1358
1359		! Open the FMS2 file-set.
1360	0	success = fms2_open_file(fileobj, filename, "read", MOM_domain%mpp_domain)
1361	0	if (.not.success) call MOM_error(FATAL, "Failed to open "//trim(filename))
1362
1363		! Find the matching case-insensitive u- and v-variable names in the file and prepare to read them.
1364		call prepare_to_read_var(fileobj, u_fieldname, "read_vector_3d: ", filename, &
1365	0	u_var, has_time_dim, timelevel, position=u_pos)
1366		call prepare_to_read_var(fileobj, v_fieldname, "read_vector_3d: ", filename, &
1367	0	v_var, has_time_dim, timelevel, position=v_pos)
1368
1369		! Read the u-data and v-data, dangerously assuming either both or neither have time dimensions.
1370		! There would already been an error message for one of the variables if they are inconsistent.
1371	0	if (present(timelevel) .and. has_time_dim) then
1372	0	call fms2_read_data(fileobj, u_var, u_data, unlim_dim_level=timelevel)
1373	0	call fms2_read_data(fileobj, v_var, v_data, unlim_dim_level=timelevel)
1374		else
1375	0	call fms2_read_data(fileobj, u_var, u_data)
1376	0	call fms2_read_data(fileobj, v_var, v_data)
1377		endif
1378
1379		! Close the file-set.
1380	0	if (check_if_open(fileobj)) call fms2_close_file(fileobj)
1381
1382	0	if (present(scale)) then ; if (scale /= 1.0) then
1383	0	call rescale_comp_data(MOM_Domain, u_data, scale)
1384	0	call rescale_comp_data(MOM_Domain, v_data, scale)
1385		endif ; endif
1386
1387	0	end subroutine read_vector_3d
1388
1389
1390		!> Find the case-sensitive name of the variable in a netCDF file with a case-insensitive name match.
1391		!! Optionally also determine whether this variable has an unlimited time dimension.
1392	0	subroutine find_varname_in_file(fileobj, fieldname, err_header, filename, var_to_read, has_time_dim, timelevel)
1393		type(FmsNetcdfFile_t), intent(inout) :: fileobj !< An FMS2 handle to an open NetCDF file
1394		character(len=*), intent(in) :: fieldname !< The variable name to seek in the file
1395		character(len=*), intent(in) :: err_header !< A descriptive prefix for error messages
1396		character(len=*), intent(in) :: filename !< The name of the file to read
1397		character(len=*), intent(out) :: var_to_read !< The variable name to read from the file
1398		logical, optional, intent(out) :: has_time_dim !< Indicates whether fieldname has a time dimension
1399		integer, optional, intent(in) :: timelevel !< A time level to read
1400
1401		! Local variables
1402		logical :: variable_found ! Is a case-insensitive version of the variable found in the netCDF file?
1403	0	character(len=256), allocatable, dimension(:) :: var_names ! The names of all the variables in the netCDF file
1404	0	character(len=256), allocatable :: dim_names(:) ! The names of a variable's dimensions
1405		integer :: nvars ! The number of variables in the file
1406		integer :: dim_unlim_size ! The current size of the unlimited (time) dimension in the file.
1407		integer :: num_var_dims ! The number of dimensions a variable has in the file.
1408		integer :: time_dim ! The position of the unlimited (time) dimension for a variable, or -1
1409		! if it has no unlimited dimension.
1410		integer :: i
1411
1412		! Open the file if necessary
1413	0	if (.not.check_if_open(fileobj)) &
1414	0	call MOM_error(FATAL, trim(err_header)//trim(filename)//" was not open in call to find_varname_in_file.")
1415
1416		! Search for the variable in the file, looking for the case-sensitive name first.
1417	0	if (variable_exists(fileobj, trim(fieldname))) then
1418	0	var_to_read = trim(fieldname)
1419		else ! Look for case-insensitive variable name matches.
1420	0	nvars = get_num_variables(fileobj)
1421	0	if (nvars < 1) call MOM_error(FATAL, "nvars is less than 1 for file "//trim(filename))
1422	0	allocate(var_names(nvars))
1423	0	call get_variable_names(fileobj, var_names)
1424
1425		! search for the variable in the file
1426	0	variable_found = .false.
1427	0	do i=1,nvars
1428	0	if (lowercase(trim(var_names(i))) == lowercase(trim(fieldname))) then
1429	0	variable_found = .true.
1430	0	var_to_read = trim(var_names(i))
1431	0	exit
1432		endif
1433		enddo
1434	0	if (.not.(variable_found)) &
1435	0	call MOM_error(FATAL, trim(err_header)//trim(fieldname)//" not found in "//trim(filename))
1436	0	deallocate(var_names)
1437		endif
1438
1439		! FMS2 can not handle a timelevel argument if the variable does not have one in the file,
1440		! so some error checking and logic are required.
1441	0	if (present(has_time_dim) .or. present(timelevel)) then
1442	0	time_dim = -1
1443
1444	0	num_var_dims = get_variable_num_dimensions(fileobj, trim(var_to_read))
1445	0	allocate(dim_names(num_var_dims)) ; dim_names(:) = ""
1446	0	call get_variable_dimension_names(fileobj, trim(var_to_read), dim_names)
1447
1448	0	do i=1,num_var_dims
1449	0	if (is_dimension_unlimited(fileobj, dim_names(i))) then
1450	0	time_dim = i
1451	0	if (present(timelevel)) then
1452	0	call get_dimension_size(fileobj, dim_names(i), dim_unlim_size)
1453	0	if ((timelevel > dim_unlim_size) .and. is_root_PE()) call MOM_error(FATAL, &
1454		trim(err_header)//"Attempting to read a time level of "//trim(var_to_read)//&
1455	0	" that exceeds the size of the time dimension in "//trim(filename))
1456		endif
1457	0	exit
1458		endif
1459		enddo
1460	0	deallocate(dim_names)
1461
1462	0	if (present(timelevel) .and. (time_dim < 0) .and. is_root_PE()) &
1463		call MOM_error(WARNING, trim(err_header)//"time level specified, but the variable "//&
1464	0	trim(var_to_read)//" does not have an unlimited dimension in "//trim(filename))
1465	0	if ((.not.present(timelevel)) .and. (time_dim > 0) .and. is_root_PE()) &
1466		call MOM_error(WARNING, trim(err_header)//"The variable "//trim(var_to_read)//&
1467	0	" has an unlimited dimension in "//trim(filename)//" but no time level is specified.")
1468	0	if (present(has_time_dim)) has_time_dim = (time_dim > 0)
1469		endif
1470
1471	0	end subroutine find_varname_in_file
1472
1473
1474		!> Find the case-insensitive name match with a variable in an open domain-decomposed file-set,
1475		!! prepare FMS2 to read this variable, and return some information needed to call fms2_read_data
1476		!! correctly for this variable and file.
1477	0	subroutine prepare_to_read_var(fileobj, fieldname, err_header, filename, var_to_read, &
1478		has_time_dim, timelevel, position)
1479		type(FmsNetcdfDomainFile_t), intent(inout) :: fileobj !< An FMS2 handle to an open domain-decomposed file
1480		character(len=*), intent(in) :: fieldname !< The variable name to seek in the file
1481		character(len=*), intent(in) :: err_header !< A descriptive prefix for error messages
1482		character(len=*), intent(in) :: filename !< The name of the file to read
1483		character(len=*), intent(out) :: var_to_read !< The variable name to read from the file
1484		logical, optional, intent(out) :: has_time_dim !< Indicates whether fieldname has a time dimension
1485		integer, optional, intent(in) :: timelevel !< A time level to read
1486		integer, optional, intent(in) :: position !< A flag indicating where this variable is discretized
1487
1488		! Local variables
1489		logical :: variable_found ! Is a case-insensitive version of the variable found in the netCDF file?
1490	0	character(len=256), allocatable, dimension(:) :: var_names ! The names of all the variables in the netCDF file
1491	0	character(len=256), allocatable :: dim_names(:) ! The names of a variable's dimensions
1492		integer :: nvars ! The number of variables in the file.
1493		integer :: dim_unlim_size ! The current size of the unlimited (time) dimension in the file.
1494		integer :: num_var_dims ! The number of dimensions a variable has in the file.
1495		integer :: time_dim ! The position of the unlimited (time) dimension for a variable, or -1
1496		! if it has no unlimited dimension.
1497		integer :: i
1498
1499		! Open the file if necessary
1500	0	if (.not.check_if_open(fileobj)) &
1501	0	call MOM_error(FATAL, trim(err_header)//trim(filename)//" was not open in call to prepare_to_read_var.")
1502
1503		! Search for the variable in the file, looking for the case-sensitive name first.
1504	0	if (variable_exists(fileobj, trim(fieldname))) then
1505	0	var_to_read = trim(fieldname)
1506		else ! Look for case-insensitive variable name matches.
1507	0	nvars = get_num_variables(fileobj)
1508	0	if (nvars < 1) call MOM_error(FATAL, "nvars is less than 1 for file "//trim(filename))
1509	0	allocate(var_names(nvars))
1510	0	call get_variable_names(fileobj, var_names)
1511
1512	0	variable_found = .false.
1513	0	do i=1,nvars
1514	0	if (lowercase(trim(var_names(i))) == lowercase(trim(fieldname))) then
1515	0	variable_found = .true.
1516	0	var_to_read = trim(var_names(i))
1517	0	exit
1518		endif
1519		enddo
1520	0	if (.not.(variable_found)) &
1521	0	call MOM_error(FATAL, trim(err_header)//trim(fieldname)//" not found in "//trim(filename))
1522	0	deallocate(var_names)
1523		endif
1524
1525		! FMS2 can not handle a timelevel argument if the variable does not have one in the file,
1526		! so some error checking and logic are required.
1527	0	if (present(has_time_dim) .or. present(timelevel)) then
1528	0	time_dim = -1
1529
1530	0	num_var_dims = get_variable_num_dimensions(fileobj, trim(var_to_read))
1531	0	allocate(dim_names(num_var_dims)) ; dim_names(:) = ""
1532	0	call get_variable_dimension_names(fileobj, trim(var_to_read), dim_names)
1533
1534	0	do i=1,num_var_dims
1535	0	if (is_dimension_unlimited(fileobj, dim_names(i))) then
1536	0	time_dim = i
1537	0	if (present(timelevel)) then
1538	0	call get_dimension_size(fileobj, dim_names(i), dim_unlim_size)
1539	0	if ((timelevel > dim_unlim_size) .and. is_root_PE()) call MOM_error(FATAL, &
1540		trim(err_header)//"Attempting to read a time level of "//trim(var_to_read)//&
1541	0	" that exceeds the size of the time dimension in "//trim(filename))
1542		endif
1543	0	exit
1544		endif
1545		enddo
1546	0	deallocate(dim_names)
1547
1548	0	if (present(timelevel) .and. (time_dim < 0) .and. is_root_PE()) &
1549		call MOM_error(WARNING, trim(err_header)//"time level specified, but the variable "//&
1550	0	trim(var_to_read)//" does not have an unlimited dimension in "//trim(filename))
1551	0	if ((.not.present(timelevel)) .and. (time_dim > 0) .and. is_root_PE()) &
1552		call MOM_error(WARNING, trim(err_header)//"The variable "//trim(var_to_read)//&
1553	0	" has an unlimited dimension in "//trim(filename)//" but no time level is specified.")
1554	0	if (present(has_time_dim)) has_time_dim = (time_dim > 0)
1555		endif
1556
1557		! Registering the variable axes essentially just specifies the discrete position of this variable.
1558	0	call MOM_register_variable_axes(fileobj, var_to_read, filename, position)
1559
1560	0	end subroutine prepare_to_read_var
1561
1562		!> register axes associated with a variable from a domain-decomposed netCDF file
1563	0	subroutine MOM_register_variable_axes(fileObj, variableName, filename, position)
1564		type(FmsNetcdfDomainFile_t), intent(inout) :: fileObj !< Handle to an open FMS2 netCDF file object
1565		character(len=*), intent(in) :: variableName !< name of the variable
1566		character(len=*), intent(in) :: filename !< The name of the file to read
1567		integer, optional, intent(in) :: position !< A flag indicating where this data is discretized
1568
1569		! Local variables
1570	0	character(len=256), allocatable, dimension(:) :: dim_names ! variable dimension names
1571	0	integer, allocatable, dimension(:) :: dimSizes ! variable dimension sizes
1572	0	logical, allocatable, dimension(:) :: is_x ! Is this a (likely domain-decomposed) x-axis
1573	0	logical, allocatable, dimension(:) :: is_y ! Is this a (likely domain-decomposed) y-axis
1574	0	logical, allocatable, dimension(:) :: is_t ! Is this a time axis or another unlimited axis
1575		integer :: ndims ! number of dimensions
1576		integer :: xPos, yPos ! Discrete positions for x and y axes. Default is CENTER
1577		integer :: i
1578
1579	0	xPos = CENTER ; yPos = CENTER
1580	0	if (present(position)) then
1581	0	if ((position == CORNER) .or. (position == EAST_FACE)) xPos = EAST_FACE
1582	0	if ((position == CORNER) .or. (position == NORTH_FACE)) yPos = NORTH_FACE
1583		endif
1584
1585		! get variable dimension names and lengths
1586	0	ndims = get_variable_num_dimensions(fileObj, trim(variableName))
1587	0	allocate(dimSizes(ndims))
1588	0	allocate(dim_names(ndims))
1589	0	allocate(is_x(ndims)) ; is_x(:) = .false.
1590	0	allocate(is_y(ndims)) ; is_y(:) = .false.
1591	0	allocate(is_t(ndims)) ; is_t(:) = .false.
1592	0	call get_variable_size(fileObj, trim(variableName), dimSizes)
1593	0	call get_variable_dimension_names(fileObj, trim(variableName), dim_names)
1594	0	call categorize_axes(fileObj, filename, ndims, dim_names, is_x, is_y, is_t)
1595
1596		! register the axes
1597	0	do i=1,ndims
1598	0	if ( .not.is_dimension_registered(fileobj, trim(dim_names(i))) ) then
1599	0	if (is_x(i)) then
1600	0	call register_axis(fileObj, trim(dim_names(i)), "x", domain_position=xPos)
1601	0	elseif (is_y(i)) then
1602	0	call register_axis(fileObj, trim(dim_names(i)), "y", domain_position=yPos)
1603		else
1604	0	call register_axis(fileObj, trim(dim_names(i)), dimSizes(i))
1605		endif
1606		endif
1607		enddo
1608
1609	0	deallocate(dimSizes, dim_names, is_x, is_y, is_t)
1610	0	end subroutine MOM_register_variable_axes
1611
1612		!> Determine whether a variable's axes are associated with x-, y- or time-dimensions. Other
1613		!! unlimited dimensions are also labeled as time axes for these purposes.
1614	0	subroutine categorize_axes(fileObj, filename, ndims, dim_names, is_x, is_y, is_t)
1615		class(FmsNetcdfFile_t), intent(in) :: fileObj !< Handle to an open FMS2 netCDF file object
1616		character(len=*), intent(in) :: filename !< The name of the file to read
1617		integer, intent(in) :: ndims !< The number of dimensions associated with a variable
1618		character(len=*), dimension(ndims), intent(in) :: dim_names !< Names of the dimensions associated with a variable
1619		logical, dimension(ndims), intent(out) :: is_x !< Indicates if each dimension a (likely decomposed) x-axis
1620		logical, dimension(ndims), intent(out) :: is_y !< Indicates if each dimension a (likely decomposed) y-axis
1621		logical, dimension(ndims), intent(out) :: is_t !< Indicates if each dimension unlimited (usually time) axis
1622
1623		! Local variables
1624		character(len=128) :: cartesian ! A flag indicating a Cartesian direction - usually a single character.
1625		character(len=512) :: dim_list ! A concatenated list of dimension names.
1626		character(len=128) :: units ! units corresponding to a specific variable dimension
1627		logical :: x_found, y_found ! Indicate whether an x- or y- dimension have been found.
1628		integer :: i
1629
1630	0	x_found = .false. ; y_found = .false.
1631	0	is_x(:) = .false. ; is_y(:) = .false.
1632	0	do i=1,ndims
1633	0	is_t(i) = is_dimension_unlimited(fileObj, trim(dim_names(i)))
1634		! First look for indicative variable attributes
1635	0	if (.not.is_t(i)) then
1636	0	if (variable_exists(fileobj, trim(dim_names(i)))) then
1637	0	cartesian = ""
1638	0	if (variable_att_exists(fileobj, trim(dim_names(i)), "cartesian_axis")) then
1639	0	call get_variable_attribute(fileobj, trim(dim_names(i)), "cartesian_axis", cartesian(1:1))
1640	0	elseif (variable_att_exists(fileobj, trim(dim_names(i)), "axis")) then
1641	0	call get_variable_attribute(fileobj, trim(dim_names(i)), "axis", cartesian(1:1))
1642		endif
1643	0	cartesian = adjustl(cartesian)
1644	0	if ((index(cartesian, "X") == 1) .or. (index(cartesian, "x") == 1)) is_x(i) = .true.
1645	0	if ((index(cartesian, "Y") == 1) .or. (index(cartesian, "y") == 1)) is_y(i) = .true.
1646	0	if ((index(cartesian, "T") == 1) .or. (index(cartesian, "t") == 1)) is_t(i) = .true.
1647		endif
1648		endif
1649	0	if (is_x(i)) x_found = .true.
1650	0	if (is_y(i)) y_found = .true.
1651		enddo
1652
1653	0	if (.not.(x_found .and. y_found)) then
1654		! Next look for hints from axis names for uncharacterized axes
1655	0	do i=1,ndims ; if (.not.(is_x(i) .or. is_y(i) .or. is_t(i))) then
1656	0	call categorize_axis_from_name(dim_names(i), is_x(i), is_y(i))
1657	0	if (is_x(i)) x_found = .true.
1658	0	if (is_y(i)) y_found = .true.
1659		endif ; enddo
1660		endif
1661
1662	0	if (.not.(x_found .and. y_found)) then
1663		! Look for hints from CF-compliant axis units for uncharacterized axes
1664	0	do i=1,ndims ; if (.not.(is_x(i) .or. is_y(i) .or. is_t(i))) then
1665	0	if (variable_exists(fileobj, trim(dim_names(i)))) then
1666	0	call get_variable_units(fileobj, trim(dim_names(i)), units)
1667	0	call categorize_axis_from_units(units, is_x(i), is_y(i))
1668		endif
1669	0	if (is_x(i)) x_found = .true.
1670	0	if (is_y(i)) y_found = .true.
1671		endif ; enddo
1672		endif
1673
1674	0	if (.not.(x_found .and. y_found) .and. ((ndims>2) .or. ((ndims==2) .and. .not.is_t(ndims)))) then
1675		! This is a case where one would expect to find x-and y-dimensions, but none have been found.
1676	0	if (is_root_pe()) then
1677	0	dim_list = trim(dim_names(1))//", "//trim(dim_names(2))
1678	0	do i=3,ndims ; dim_list = trim(dim_list)//", "//trim(dim_names(i)) ; enddo
1679		call MOM_error(WARNING, "categorize_axes: Failed to identify x- and y- axes in the axis list ("//&
1680	0	trim(dim_list)//") of a variable being read from "//trim(filename))
1681		endif
1682		endif
1683
1684	0	end subroutine categorize_axes
1685
1686		!> Determine whether an axis is associated with the x- or y-directions based on a comparison of
1687		!! its units with CF-compliant variants of latitude or longitude units.
1688	0	subroutine categorize_axis_from_units(unit_string, is_x, is_y)
1689		character(len=*), intent(in) :: unit_string !< string of units
1690		logical, intent(out) :: is_x !< Indicates if the axis units are associated with an x-direction axis
1691		logical, intent(out) :: is_y !< Indicates if the axis units are associated with an y-direction axis
1692
1693	0	is_x = .false. ; is_y = .false.
1694	0	select case (lowercase(trim(unit_string)))
1695	0	case ("degrees_north"); is_y = .true.
1696	0	case ("degree_north") ; is_y = .true.
1697	0	case ("degrees_n") ; is_y = .true.
1698	0	case ("degree_n") ; is_y = .true.
1699	0	case ("degreen") ; is_y = .true.
1700	0	case ("degreesn") ; is_y = .true.
1701	0	case ("degrees_east") ; is_x = .true.
1702	0	case ("degree_east") ; is_x = .true.
1703	0	case ("degreese") ; is_x = .true.
1704	0	case ("degreee") ; is_x = .true.
1705	0	case ("degree_e") ; is_x = .true.
1706	0	case ("degrees_e") ; is_x = .true.
1707	0	case default ; is_x = .false. ; is_y = .false.
1708		end select
1709
1710	0	end subroutine categorize_axis_from_units
1711
1712		!> Tries to determine whether the axis name is commonly associated with an x- or y- axis. This
1713		!! approach is fragile and unreliable, but it a backup to reading a CARTESIAN file attribute.
1714	0	subroutine categorize_axis_from_name(dimname, is_x, is_y)
1715		character(len=*), intent(in) :: dimname !< A dimension name
1716		logical, intent(out) :: is_x !< Indicates if the axis name is associated with an x-direction axis
1717		logical, intent(out) :: is_y !< Indicates if the axis name is associated with an y-direction axis
1718
1719	0	is_x = .false. ; is_y = .false.
1720	0	select case(trim(lowercase(dimname)))
1721	0	case ("grid_x_t") ; is_x = .true.
1722	0	case ("nx") ; is_x = .true.
1723	0	case ("nxp") ; is_x = .true.
1724	0	case ("longitude") ; is_x = .true.
1725	0	case ("long") ; is_x = .true.
1726	0	case ("lon") ; is_x = .true.
1727	0	case ("lonh") ; is_x = .true.
1728	0	case ("lonq") ; is_x = .true.
1729	0	case ("xh") ; is_x = .true.
1730	0	case ("xq") ; is_x = .true.
1731	0	case ("i") ; is_x = .true.
1732
1733	0	case ("grid_y_t") ; is_y = .true.
1734	0	case ("ny") ; is_y = .true.
1735	0	case ("nyp") ; is_y = .true.
1736	0	case ("latitude") ; is_y = .true.
1737	0	case ("lat") ; is_y = .true.
1738	0	case ("lath") ; is_y = .true.
1739	0	case ("latq") ; is_y = .true.
1740	0	case ("yh") ; is_y = .true.
1741	0	case ("yq") ; is_y = .true.
1742	0	case ("j") ; is_y = .true.
1743
1744	0	case default ; is_x = .false. ; is_y = .false.
1745		end select
1746
1747	0	end subroutine categorize_axis_from_name
1748
1749
1750		!> Write a 4d field to an output file.
1751	0	subroutine write_field_4d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, fill_value)
1752		type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing
1753		type(fieldtype), intent(in) :: field_md !< Field type with metadata
1754		type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition
1755		real, dimension(:,:,:,:), intent(inout) :: field !< Field to write
1756		real, optional, intent(in) :: tstamp !< Model time of this field
1757		integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1)
1758		real, optional, intent(in) :: fill_value !< Missing data fill value
1759
1760
1761		! Local variables
1762		integer :: time_index
1763
1764	0	if (present(tstamp)) then
1765	0	time_index = write_time_if_later(IO_handle, tstamp)
1766	0	call write_data(IO_handle%fileobj, trim(field_md%name), field, unlim_dim_level=time_index)
1767		else
1768	0	call write_data(IO_handle%fileobj, trim(field_md%name), field)
1769		endif
1770	0	end subroutine write_field_4d
1771
1772		!> Write a 3d field to an output file.
1773	31	subroutine write_field_3d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, fill_value)
1774		type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing
1775		type(fieldtype), intent(in) :: field_md !< Field type with metadata
1776		type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition
1777		real, dimension(:,:,:), intent(inout) :: field !< Field to write
1778		real, optional, intent(in) :: tstamp !< Model time of this field
1779		integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1)
1780		real, optional, intent(in) :: fill_value !< Missing data fill value
1781
1782		! Local variables
1783		integer :: time_index
1784
1785	31	if (present(tstamp)) then
1786	31	time_index = write_time_if_later(IO_handle, tstamp)
1787	31	call write_data(IO_handle%fileobj, trim(field_md%name), field, unlim_dim_level=time_index)
1788		else
1789	0	call write_data(IO_handle%fileobj, trim(field_md%name), field)
1790		endif
1791	31	end subroutine write_field_3d
1792
1793		!> Write a 2d field to an output file.
1794	28	subroutine write_field_2d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, fill_value)
1795		type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing
1796		type(fieldtype), intent(in) :: field_md !< Field type with metadata
1797		type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition
1798		real, dimension(:,:), intent(inout) :: field !< Field to write
1799		real, optional, intent(in) :: tstamp !< Model time of this field
1800		integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1)
1801		real, optional, intent(in) :: fill_value !< Missing data fill value
1802
1803		! Local variables
1804		integer :: time_index
1805
1806	28	if (present(tstamp)) then
1807	28	time_index = write_time_if_later(IO_handle, tstamp)
1808	28	call write_data(IO_handle%fileobj, trim(field_md%name), field, unlim_dim_level=time_index)
1809		else
1810	0	call write_data(IO_handle%fileobj, trim(field_md%name), field)
1811		endif
1812	28	end subroutine write_field_2d
1813
1814		!> Write a 1d field to an output file.
1815	0	subroutine write_field_1d(IO_handle, field_md, field, tstamp)
1816		type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing
1817		type(fieldtype), intent(in) :: field_md !< Field type with metadata
1818		real, dimension(:), intent(in) :: field !< Field to write
1819		real, optional, intent(in) :: tstamp !< Model time of this field
1820
1821		! Local variables
1822		integer :: time_index
1823
1824	0	if (present(tstamp)) then
1825	0	time_index = write_time_if_later(IO_handle, tstamp)
1826	0	call write_data(IO_handle%fileobj, trim(field_md%name), field, unlim_dim_level=time_index)
1827		else
1828	0	call write_data(IO_handle%fileobj, trim(field_md%name), field)
1829		endif
1830	0	end subroutine write_field_1d
1831
1832		!> Write a 0d field to an output file.
1833	4	subroutine write_field_0d(IO_handle, field_md, field, tstamp)
1834		type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing
1835		type(fieldtype), intent(in) :: field_md !< Field type with metadata
1836		real, intent(in) :: field !< Field to write
1837		real, optional, intent(in) :: tstamp !< Model time of this field
1838
1839		! Local variables
1840		integer :: time_index
1841
1842	4	if (present(tstamp)) then
1843	4	time_index = write_time_if_later(IO_handle, tstamp)
1844	4	call write_data(IO_handle%fileobj, trim(field_md%name), field, unlim_dim_level=time_index)
1845		else
1846	0	call write_data(IO_handle%fileobj, trim(field_md%name), field)
1847		endif
1848	4	end subroutine write_field_0d
1849
1850		!> Returns the integer time index for a write in this file, also writing the time variable to
1851		!! the file if this time is later than what is already in the file.
1852	63	integer function write_time_if_later(IO_handle, field_time)
1853		type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing
1854		real, intent(in) :: field_time !< Model time of this field
1855
1856		! Local variables
1857		character(len=256) :: dim_unlim_name ! name of the unlimited dimension in the file
1858
1859	63	if ((field_time > IO_handle%file_time) .or. (IO_handle%num_times == 0)) then
1860	2	IO_handle%file_time = field_time
1861	2	IO_handle%num_times = IO_handle%num_times + 1
1862	2	dim_unlim_name = find_unlimited_dimension_name(IO_handle%fileobj)
1863	2	if (len_trim(dim_unlim_name) > 0) &
1864		call write_data(IO_handle%fileobj, trim(dim_unlim_name), [field_time], &
1865	6	corner=[IO_handle%num_times], edge_lengths=[1])
1866		endif
1867
1868	63	write_time_if_later = IO_handle%num_times
1869	63	end function write_time_if_later
1870
1871		!> Write the data for an axis
1872	12	subroutine MOM_write_axis(IO_handle, axis)
1873		type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing
1874		type(axistype), intent(in) :: axis !< An axis type variable with information to write
1875
1876		integer :: is, ie
1877
1878	12	if (axis%domain_decomposed) then
1879		! FMS2 does not domain-decompose 1d arrays, so we explicitly slice it
1880	8	call get_global_io_domain_indices(IO_handle%fileobj, trim(axis%name), is, ie)
1881	8	call write_data(IO_handle%fileobj, trim(axis%name), axis%ax_data(is:ie))
1882		else
1883	4	call write_data(IO_handle%fileobj, trim(axis%name), axis%ax_data)
1884		endif
1885	12	end subroutine MOM_write_axis
1886
1887		!> Store information about an axis in a previously defined axistype and write this
1888		!! information to the file indicated by unit.
1889	14	subroutine write_metadata_axis(IO_handle, axis, name, units, longname, cartesian, sense, domain, data, &
1890		edge_axis, calendar)
1891		type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing
1892		type(axistype), intent(inout) :: axis !< The axistype where this information is stored.
1893		character(len=*), intent(in) :: name !< The name in the file of this axis
1894		character(len=*), intent(in) :: units !< The units of this axis
1895		character(len=*), intent(in) :: longname !< The long description of this axis
1896		character(len=*), optional, intent(in) :: cartesian !< A variable indicating which direction
1897		!! this axis corresponds with. Valid values
1898		!! include 'X', 'Y', 'Z', 'T', and 'N' for none.
1899		integer, optional, intent(in) :: sense !< This is 1 for axes whose values increase upward, or
1900		!! -1 if they increase downward.
1901		type(domain1D), optional, intent(in) :: domain !< The domain decomposion for this axis
1902		real, dimension(:), optional, intent(in) :: data !< The coordinate values of the points on this axis
1903		logical, optional, intent(in) :: edge_axis !< If true, this axis marks an edge of the tracer cells
1904		character(len=*), optional, intent(in) :: calendar !< The name of the calendar used with a time axis
1905
1906	14	character(len=:), allocatable :: cart ! A left-adjusted and trimmed copy of cartesian
1907		logical :: is_x, is_y, is_t ! If true, this is a domain-decomposed axis in one of the directions.
1908		integer :: position ! A flag indicating the axis staggering position.
1909		integer :: i, isc, iec, global_size
1910
1911	14	if (is_dimension_registered(IO_handle%fileobj, trim(name))) then
1912		call MOM_error(FATAL, "write_metadata_axis was called more than once for axis "//trim(name)//&
1913	0	" in file "//trim(IO_handle%filename))
1914	0	return
1915		endif
1916
1917	14	axis%name = trim(name)
1918	14	if (present(data) .and. allocated(axis%ax_data)) call MOM_error(FATAL, &
1919		"Data is already allocated in a call to write_metadata_axis for axis "//&
1920	0	trim(name)//" in file "//trim(IO_handle%filename))
1921
1922	14	is_x = .false. ; is_y = .false. ; is_t = .false.
1923	14	position = CENTER
1924	14	if (present(cartesian)) then
1925	14	cart = trim(adjustl(cartesian))
1926	14	if ((index(cart, "X") == 1) .or. (index(cart, "x") == 1)) is_x = .true.
1927	14	if ((index(cart, "Y") == 1) .or. (index(cart, "y") == 1)) is_y = .true.
1928	14	if ((index(cart, "T") == 1) .or. (index(cart, "t") == 1)) is_t = .true.
1929		endif
1930
1931		! For now, we assume that all horizontal axes are domain-decomposed.
1932	14	if (is_x .or. is_y) &
1933	8	axis%domain_decomposed = .true.
1934
1935	14	if (is_x) then
1936	4	if (present(edge_axis)) then ; if (edge_axis) position = EAST_FACE ; endif
1937	4	call register_axis(IO_handle%fileobj, trim(name), 'x', domain_position=position)
1938	10	elseif (is_y) then
1939	4	if (present(edge_axis)) then ; if (edge_axis) position = NORTH_FACE ; endif
1940	4	call register_axis(IO_handle%fileobj, trim(name), 'y', domain_position=position)
1941	6	elseif (is_t .and. .not.present(data)) then
1942		! This is the unlimited (time) dimension.
1943	2	call register_axis(IO_handle%fileobj, trim(name), unlimited)
1944		else
1945	4	if (.not.present(data)) call MOM_error(FATAL,"MOM_io:register_diagnostic_axis: "//&
1946	0	"An axis_length argument is required to register the axis "//trim(name))
1947	4	call register_axis(IO_handle%fileobj, trim(name), size(data))
1948		endif
1949
1950	14	if (present(data)) then
1951		! With FMS2, the data for the axis labels has to match the computational domain on this PE.
1952	12	if (present(domain)) then
1953		! The commented-out code on the next ~11 lines runs but there is missing data in the output file
1954		! call mpp_get_compute_domain(domain, isc, iec)
1955		! call mpp_get_global_domain(domain, size=global_size)
1956		! if (size(data) == global_size) then
1957		! allocate(axis%ax_data(iec+1-isc)) ; axis%ax_data(:) = data(isc:iec)
1958		! ! A simpler set of labels: do i=1,iec-isc ; axis%ax_data(i) = real(isc + i) - 1.0 ; enddo
1959		! elseif (size(data) == global_size+1) then
1960		! ! This is an edge axis. Note the effective SW indexing convention here.
1961		! allocate(axis%ax_data(iec+2-isc)) ; axis%ax_data(:) = data(isc:iec+1)
1962		! ! A simpler set of labels: do i=1,iec+1-isc ; axis%ax_data(i) = real(isc + i) - 1.5 ; enddo
1963		! else
1964		! call MOM_error(FATAL, "Unexpected size of data for "//trim(name)//" in write_metadata_axis.")
1965		! endif
1966
1967		! This works for a simple 1x1 IO layout, but gives errors for nontrivial IO layouts
1968	732	allocate(axis%ax_data(size(data))) ; axis%ax_data(:) = data(:)
1969
1970		else ! Store the entire array of axis labels.
1971	306	allocate(axis%ax_data(size(data))) ; axis%ax_data(:) = data(:)
1972		endif
1973		endif
1974
1975
1976		! Now create the variable that describes this axis.
1977	28	call register_field(IO_handle%fileobj, trim(name), "double", dimensions=(/name/))
1978	14	if (len_trim(longname) > 0) &
1979		call register_variable_attribute(IO_handle%fileobj, trim(name), 'long_name', &
1980	14	trim(longname), len_trim(longname))
1981	14	if (len_trim(units) > 0) &
1982		call register_variable_attribute(IO_handle%fileobj, trim(name), 'units', &
1983	14	trim(units), len_trim(units))
1984	14	if (present(cartesian)) &
1985		call register_variable_attribute(IO_handle%fileobj, trim(name), 'cartesian_axis', &
1986	14	trim(cartesian), len_trim(cartesian))
1987	14	if (present(sense)) &
1988	4	call register_variable_attribute(IO_handle%fileobj, trim(name), 'sense', sense)
1989	28	end subroutine write_metadata_axis
1990
1991		!> Store information about an output variable in a previously defined fieldtype and write this
1992		!! information to the file indicated by unit.
1993	63	subroutine write_metadata_field(IO_handle, field, axes, name, units, longname, &
1994	63	pack, standard_name, checksum)
1995		type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing
1996		type(fieldtype), intent(inout) :: field !< The fieldtype where this information is stored
1997		type(axistype), dimension(:), intent(in) :: axes !< Handles for the axis used for this variable
1998		character(len=*), intent(in) :: name !< The name in the file of this variable
1999		character(len=*), intent(in) :: units !< The units of this variable
2000		character(len=*), intent(in) :: longname !< The long description of this variable
2001		integer, optional, intent(in) :: pack !< A precision reduction factor with which the
2002		!! variable. The default, 1, has no reduction,
2003		!! but 2 is not uncommon.
2004		character(len=*), optional, intent(in) :: standard_name !< The standard (e.g., CMOR) name for this variable
2005		integer(kind=int64), dimension(:), &
2006		optional, intent(in) :: checksum !< Checksum values that can be used to verify reads.
2007
2008		! Local variables
2009	126	character(len=256), dimension(size(axes)) :: dim_names ! The names of the dimensions
2010		character(len=16) :: prec_string ! A string specifying the precision with which to save this variable
2011		character(len=64) :: checksum_string ! checksum character array created from checksum argument
2012		integer :: i, ndims
2013
2014	63	ndims = size(axes)
2015	275	do i=1,ndims ; dim_names(i) = trim(axes(i)%name) ; enddo
2016	63	prec_string = "double" ; if (present(pack)) then ; if (pack > 1) prec_string = "float" ; endif
2017	63	call register_field(IO_handle%fileobj, trim(name), trim(prec_string), dimensions=dim_names)
2018	63	if (len_trim(longname) > 0) &
2019		call register_variable_attribute(IO_handle%fileobj, trim(name), 'long_name', &
2020	63	trim(longname), len_trim(longname))
2021	63	if (len_trim(units) > 0) &
2022		call register_variable_attribute(IO_handle%fileobj, trim(name), 'units', &
2023	63	trim(units), len_trim(units))
2024	63	if (present(standard_name)) &
2025		call register_variable_attribute(IO_handle%fileobj, trim(name), 'standard_name', &
2026	0	trim(standard_name), len_trim(standard_name))
2027	63	if (present(checksum)) then
2028	63	write (checksum_string,'(Z16)') checksum(1) ! Z16 is the hexadecimal format code
2029		call register_variable_attribute(IO_handle%fileobj, trim(name), "checksum", &
2030	63	trim(checksum_string), len_trim(checksum_string))
2031		endif
2032
2033		! Store information in the field-type, regardless of which interfaces are used.
2034	63	field%name = trim(name)
2035	63	field%longname = trim(longname)
2036	63	field%units = trim(units)
2037	63	field%chksum_read = -1
2038	63	field%valid_chksum = .false.
2039
2040	63	end subroutine write_metadata_field
2041
2042		!> Write a global text attribute to a file.
2043	0	subroutine write_metadata_global(IO_handle, name, attribute)
2044		type(file_type), intent(in) :: IO_handle !< Handle for a file that is open for writing
2045		character(len=*), intent(in) :: name !< The name in the file of this global attribute
2046		character(len=*), intent(in) :: attribute !< The value of this attribute
2047
2048	0	call register_global_attribute(IO_handle%fileobj, name, attribute, len_trim(attribute))
2049	0	end subroutine write_metadata_global
2050
2051		!> Return unlimited dimension name in file, or empty string if none exists.
2052	2	function find_unlimited_dimension_name(fileobj) result(label)
2053		type(FmsNetcdfDomainFile_t), intent(in) :: fileobj
2054		!< File handle
2055		character(len=:), allocatable :: label
2056		!< Unlimited dimension name, or empty string if none exists
2057
2058		integer :: ndims
2059		!< Number of dimensions
2060	2	character(len=256), allocatable :: dim_names(:)
2061		!< File handle dimension names
2062		integer :: i
2063		!< Loop index
2064
2065	2	ndims = get_num_dimensions(fileobj)
2066	2	allocate(dim_names(ndims))
2067	2	call get_dimension_names(fileobj, dim_names)
2068
2069	14	do i = 1, ndims
2070	14	if (is_dimension_unlimited(fileobj, dim_names(i))) then
2071	2	label = trim(dim_names(i))
2072	2	exit
2073		endif
2074		enddo
2075	2	deallocate(dim_names)
2076
2077	2	if (.not. allocated(label)) &
2078	0	label = ''
2079	2	end function find_unlimited_dimension_name
2080
2081		! NOTE: `lowercase is duplicated from `src/framework/MOM_string_functions.F90`
2082		! in order to avoid any dependency of the infra on the framework.
2083
2084		!> Return a string in which all uppercase letters have been replaced by
2085		!! their lowercase counterparts.
2086	0	function lowercase(input_string)
2087		character(len=*), intent(in) :: input_string !< The string to modify
2088		character(len=len(input_string)) :: lowercase !< The modified output string
2089		! This function returns a string in which all uppercase letters have been
2090		! replaced by their lowercase counterparts. It is loosely based on the
2091		! lowercase function in mpp_util.F90.
2092		integer, parameter :: co=iachar('a')-iachar('A') ! case offset
2093		integer :: k
2094
2095	0	lowercase = input_string
2096	0	do k=1, len_trim(input_string)
2097	0	if (lowercase(k:k) >= 'A' .and. lowercase(k:k) <= 'Z') &
2098	0	lowercase(k:k) = achar(ichar(lowercase(k:k))+co)
2099		enddo
2100	0	end function lowercase
2101
2102	0	end module MOM_io_infra