src/framework/MOM_io.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 I/O framework code
6		module MOM_io
7
8		use MOM_array_transform, only : allocate_rotated_array, rotate_array
9		use MOM_array_transform, only : rotate_array_pair, rotate_vector
10		use MOM_domains, only : MOM_domain_type, domain1D, broadcast, get_domain_components
11		use MOM_domains, only : rescale_comp_data, num_PEs, AGRID, BGRID_NE, CGRID_NE
12		use MOM_dyn_horgrid, only : dyn_horgrid_type
13		use MOM_ensemble_manager, only : get_ensemble_id
14		use MOM_error_handler, only : MOM_error, NOTE, FATAL, WARNING, is_root_PE
15		use MOM_file_parser, only : log_version, param_file_type
16		use MOM_grid, only : ocean_grid_type
17		use MOM_io_infra, only : read_field, read_vector
18		use MOM_io_infra, only : read_data => read_field ! Deprecated
19		use MOM_io_infra, only : read_field_chksum
20		use MOM_io_infra, only : file_exists
21		use MOM_io_infra, only : open_ASCII_file, close_file, file_is_open
22		use MOM_io_infra, only : get_field_size, field_exists, get_field_atts
23		use MOM_io_infra, only : get_axis_data, get_filename_suffix
24		use MOM_io_infra, only : write_version
25		use MOM_io_infra, only : MOM_namelist_file, check_namelist_error, io_infra_init, io_infra_end
26		use MOM_io_infra, only : APPEND_FILE, ASCII_FILE, MULTIPLE, NETCDF_FILE, OVERWRITE_FILE
27		use MOM_io_infra, only : READONLY_FILE, SINGLE_FILE, WRITEONLY_FILE
28		use MOM_io_infra, only : CENTER, CORNER, NORTH_FACE, EAST_FACE
29		use MOM_io_file, only : MOM_file, MOM_infra_file, MOM_netcdf_file
30		use MOM_io_file, only : MOM_axis, MOM_field
31		use MOM_string_functions, only : lowercase, slasher
32		use MOM_verticalGrid, only : verticalGrid_type
33
34		use iso_fortran_env, only : int32, int64, stdout_iso=>output_unit, stderr_iso=>error_unit
35		use netcdf, only : NF90_open, NF90_inq_varid, NF90_inq_varids, NF90_inquire, NF90_close
36		use netcdf, only : NF90_inquire_variable, NF90_get_var, NF90_get_att, NF90_inquire_attribute
37		use netcdf, only : NF90_strerror, NF90_inquire_dimension
38		use netcdf, only : NF90_NOWRITE, NF90_NOERR, NF90_GLOBAL, NF90_ENOTATT, NF90_CHAR
39
40		! The following are not used in MOM6, but may be used by externals (e.g. SIS2).
41		use MOM_io_infra, only : axistype ! still used but soon to be nuked
42		use MOM_io_infra, only : fieldtype
43		use MOM_io_infra, only : file_type
44		use MOM_io_infra, only : get_file_info
45		use MOM_io_infra, only : get_file_fields
46		use MOM_io_infra, only : get_file_times
47		use MOM_io_infra, only : open_file
48		use MOM_io_infra, only : write_field
49
50		implicit none ; private
51
52		! These interfaces are actually implemented in this file.
53		public :: create_MOM_file, reopen_MOM_file, cmor_long_std, ensembler, MOM_io_init
54		public :: MOM_field
55		public :: MOM_write_field, var_desc, modify_vardesc, query_vardesc, position_from_horgrid
56		public :: open_namelist_file, check_namelist_error, check_nml_error
57		public :: get_var_sizes, verify_variable_units, num_timelevels, read_variable, read_attribute
58		public :: open_file_to_read, close_file_to_read
59		! The following are simple pass throughs of routines from MOM_io_infra or other modules.
60		public :: file_exists, open_ASCII_file, close_file
61		public :: MOM_file, MOM_infra_file, MOM_netcdf_file
62		public :: field_exists, get_filename_appendix
63		public :: fieldtype, field_size, get_field_atts
64		public :: axistype, get_axis_data
65		public :: MOM_read_data, MOM_read_vector, read_field_chksum
66		public :: read_netCDF_data
67		public :: slasher, write_version_number
68		public :: io_infra_init, io_infra_end
69		public :: stdout_if_root
70		public :: get_var_axes_info
71		public :: get_axis_info
72		! This is used to set up information descibing non-domain-decomposed axes.
73		public :: axis_info, set_axis_info, delete_axis_info
74		! This is used to set up global file attributes
75		public :: attribute_info, set_attribute_info, delete_attribute_info
76		! This API is here just to support potential use by non-FMS drivers, and should not persist.
77		public :: read_data
78		!> These encoding constants are used to indicate the file format
79		public :: ASCII_FILE, NETCDF_FILE
80		!> These encoding constants are used to indicate whether the file is domain decomposed
81		public :: MULTIPLE, SINGLE_FILE
82		!> These encoding constants are used to indicate the access mode for a file
83		public :: APPEND_FILE, OVERWRITE_FILE, READONLY_FILE, WRITEONLY_FILE
84		!> These encoding constants are used to indicate the discretization position of a variable
85		public :: CENTER, CORNER, NORTH_FACE, EAST_FACE
86
87		! The following are not used in MOM6, but may be used by externals (e.g. SIS2).
88		public :: create_file
89		public :: reopen_file
90		public :: file_type
91		public :: open_file
92		public :: get_file_info
93		public :: get_file_fields
94		public :: get_file_times
95
96		!> Read a field from file using the infrastructure I/O.
97		interface MOM_read_data
98		module procedure MOM_read_data_0d
99		module procedure MOM_read_data_0d_int
100		module procedure MOM_read_data_1d
101		module procedure MOM_read_data_1d_int
102		module procedure MOM_read_data_2d
103		module procedure MOM_read_data_2d_region
104		module procedure MOM_read_data_3d
105		module procedure MOM_read_data_3d_region
106		module procedure MOM_read_data_4d
107		end interface MOM_read_data
108
109		!> Read a vector from file using the infrastructure I/O.
110		interface MOM_read_vector
111		module procedure MOM_read_vector_2d
112		module procedure MOM_read_vector_3d
113		end interface MOM_read_vector
114
115		!> Read a field using native netCDF I/O
116		!!
117		!! This function is primarily used for unstructured data which may contain
118		!! content that cannot be parsed by infrastructure I/O.
119		interface read_netCDF_data
120		! NOTE: Only 2D I/O is currently used; this should be expanded as needed.
121		module procedure read_netCDF_data_2d
122		end interface read_netCDF_data
123
124		!> Write a registered field to an output file, potentially with rotation
125		interface MOM_write_field
126		module procedure MOM_write_field_legacy_4d
127		module procedure MOM_write_field_legacy_3d
128		module procedure MOM_write_field_legacy_2d
129		module procedure MOM_write_field_legacy_1d
130		module procedure MOM_write_field_legacy_0d
131		module procedure MOM_write_field_4d
132		module procedure MOM_write_field_3d
133		module procedure MOM_write_field_2d
134		module procedure MOM_write_field_1d
135		module procedure MOM_write_field_0d
136		end interface MOM_write_field
137
138		!> Read an entire named variable from a named netCDF file using netCDF calls directly, rather
139		!! than any infrastructure routines and broadcast it from the root PE to the other PEs.
140		interface read_variable
141		module procedure read_variable_0d, read_variable_0d_int
142		module procedure read_variable_1d, read_variable_1d_int
143		module procedure read_variable_2d, read_variable_3d
144		end interface read_variable
145
146		!> Read a global or variable attribute from a named netCDF file using netCDF calls
147		!! directly, in some cases reading from the root PE before broadcasting to the other PEs.
148		interface read_attribute
149		module procedure read_attribute_str, read_attribute_real
150		module procedure read_attribute_int32, read_attribute_int64
151		end interface read_attribute
152
153		!> Type that stores information that can be used to create a non-decomposed axis.
154		type :: axis_info
155		character(len=32) :: name = "" !< The name of this axis for use in files
156		character(len=256) :: longname = "" !< A longer name describing this axis
157		character(len=48) :: units = "" !< The units of the axis labels
158		character(len=8) :: cartesian = "N" !< A variable indicating which direction
159		!! this axis corresponds with. Valid values
160		!! include 'X', 'Y', 'Z', 'T', and 'N' for none.
161		integer :: sense = 0 !< This is 1 for axes whose values increase upward, or -1
162		!! if they increase downward. The default, 0, is ignored.
163		integer :: ax_size = 0 !< The number of elements in this axis
164		real, allocatable, dimension(:) :: ax_data !< The values of the data on the axis [arbitrary]
165		end type axis_info
166
167		!> Type for describing a 3-d variable for output
168		type, public :: vardesc
169		character(len=64) :: name !< Variable name in a NetCDF file
170		character(len=48) :: units !< Physical dimensions of the variable
171		character(len=240) :: longname !< Long name of the variable
172		character(len=8) :: hor_grid !< Horizontal grid: u, v, h, q, Cu, Cv, T, Bu, or 1
173		character(len=8) :: z_grid !< Vertical grid: L, i, or 1
174		character(len=8) :: t_grid !< Time description: s, p, or 1
175		character(len=64) :: cmor_field_name !< CMOR name
176		character(len=64) :: cmor_units !< CMOR physical dimensions of the variable
177		character(len=240) :: cmor_longname !< CMOR long name of the variable
178		real :: conversion !< for unit conversions, such as needed to convert
179		!! from intensive to extensive [various] or [a A-1 ~> 1]
180		!! to undo internal dimensional rescaling
181		character(len=32) :: dim_names(5) !< The names in the file of the axes for this variable
182		integer :: position = -1 !< An integer encoding the horizontal position, it may
183		!! CENTER, CORNER, EAST_FACE, NORTH_FACE, or 0.
184		type(axis_info) :: extra_axes(5) !< dimensions other than space-time
185		end type vardesc
186
187		!> Type that stores for a global file attribute
188		type :: attribute_info ; private
189		character(len=:), allocatable :: name !< The name of this attribute
190		character(len=:), allocatable :: att_val !< The values of this attribute
191		end type attribute_info
192
193		integer, public :: stdout = stdout_iso !< standard output unit
194		integer, public :: stderr = stderr_iso !< standard output unit
195
196		! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
197		! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
198		! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
199		! vary with the Boussinesq approximation, the Boussinesq variant is given first.
200		! The functions in this module work with variables with arbitrary units, in which case the
201		! arbitrary rescaled units are indicated with [A ~> a], while the unscaled units are just [a].
202
203		contains
204
205		!> `create_MOM_file` wrapper for the legacy file handle, `file_type`.
206		!! NOTE: This function may be removed in a future release.
207	0	subroutine create_file(IO_handle, filename, vars, novars, fields, threading, &
208	0	timeunit, G, dG, GV, checksums, extra_axes, global_atts)
209		type(file_type), intent(inout) :: IO_handle
210		!< Handle for a files or fileset that is to be opened or reopened for
211		!! writing
212		character(len=*), intent(in) :: filename
213		!< full path to the file to create
214		type(vardesc), intent(in) :: vars(:)
215		!< structures describing fields written to filename
216		integer, intent(in) :: novars
217		!< number of fields written to filename
218		type(fieldtype), intent(inout) :: fields(:)
219		!< array of fieldtypes for each variable
220		integer, optional, intent(in) :: threading
221		!< SINGLE_FILE or MULTIPLE
222		real, optional, intent(in) :: timeunit
223		!< length of the units for time [s]. The default value is 86400.0, for 1
224		!! day.
225		type(ocean_grid_type), optional, intent(in) :: G
226		!< ocean horizontal grid structure; G or dG is required if the new file
227		!! uses any horizontal grid axes.
228		type(dyn_horgrid_type), optional, intent(in) :: dG
229		!< dynamic horizontal grid structure; G or dG is required if the new file
230		!! uses any horizontal grid axes.
231		type(verticalGrid_type), optional, intent(in) :: GV
232		!< ocean vertical grid structure, which is ! required if the new file uses
233		!! any vertical grid axes.
234		integer(kind=int64), optional, intent(in) :: checksums(:,:)
235		!< checksums of vars
236		type(axis_info), optional, intent(in) :: extra_axes(:)
237		!< Types with information about some axes that might be used in this file
238		type(attribute_info), optional, intent(in) :: global_atts(:)
239		!< Global attributes to write to this file
240
241	0	type(MOM_infra_file) :: new_file
242	0	type(MOM_field) :: new_fields(novars)
243
244	0	new_file%handle_infra = IO_handle
245
246		call create_MOM_file(new_file, filename, vars, novars, new_fields, &
247		threading=threading, timeunit=timeunit, G=G, dG=dG, GV=GV, &
248	0	checksums=checksums, extra_axes=extra_axes, global_atts=global_atts)
249
250	0	IO_handle = new_file%handle_infra
251	0	call new_file%get_file_fieldtypes(fields(:novars))
252	0	end subroutine create_file
253
254
255		!! Create a new netCDF file and register the MOM_fields to be written.
256	8	subroutine create_MOM_file(IO_handle, filename, vars, novars, fields, &
257	12	threading, timeunit, G, dG, GV, checksums, extra_axes, global_atts)
258		class(MOM_file), intent(inout) :: IO_handle !< Handle for a files or fileset that is to be
259		!! opened or reopened for writing
260		character(len=*), intent(in) :: filename !< full path to the file to create
261		type(vardesc), intent(in) :: vars(:) !< structures describing fields written to filename
262		integer, intent(in) :: novars !< number of fields written to filename
263		type(MOM_field), intent(inout) :: fields(:) !< array of fieldtypes for each variable
264		integer, optional, intent(in) :: threading !< SINGLE_FILE or MULTIPLE
265		real, optional, intent(in) :: timeunit !< length of the units for time [s]. The
266		!! default value is 86400.0, for 1 day.
267		type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG
268		!! is required if the new file uses any
269		!! horizontal grid axes.
270		type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG
271		!! is required if the new file uses any
272		!! horizontal grid axes.
273		type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is
274		!! required if the new file uses any
275		!! vertical grid axes.
276		integer(kind=int64), optional, intent(in) :: checksums(:,:) !< checksums of vars
277		type(axis_info), dimension(:), &
278		optional, intent(in) :: extra_axes !< Types with information about
279		!! some axes that might be used in this file
280		type(attribute_info), optional, intent(in) :: global_atts(:) !< Global attributes to
281		!! write to this file
282
283		logical :: use_lath, use_lonh, use_latq, use_lonq, use_time
284		logical :: use_layer, use_int, use_periodic
285		logical :: one_file, domain_set, dim_found
286	4	logical, dimension(:), allocatable :: use_extra_axis
287	4	type(MOM_axis) :: axis_lath, axis_latq, axis_lonh, axis_lonq
288	4	type(MOM_axis) :: axis_layer, axis_int, axis_time, axis_periodic
289	4	type(MOM_axis), dimension(:), allocatable :: more_axes ! Axes generated from extra_axes
290	68	type(MOM_axis) :: axes(5) ! The axes of a variable
291		type(MOM_domain_type), pointer :: Domain => NULL()
292		type(domain1d) :: x_domain, y_domain
293		integer :: position, numaxes, pack, thread, k, n, m
294		integer :: num_extra_dims ! The number of extra possible dimensions from extra_axes
295		integer :: isg, ieg, jsg, jeg, IsgB, IegB, JsgB, JegB
296		integer :: var_periods, num_periods=0
297	4	real, dimension(:), allocatable :: axis_val ! Axis label values [various]
298		real, pointer, dimension(:) :: &
299		gridLatT => NULL(), & ! The latitude of T or B points for the purpose of labeling
300		gridLatB => NULL(), & ! the output axes, often in units of [degrees_N] or [km] or [m].
301		gridLonT => NULL(), & ! The longitude of T or B points for the purpose of labeling
302		gridLonB => NULL() ! the output axes, often in units of [degrees_E] or [km] or [m].
303		character(len=40) :: time_units, x_axis_units, y_axis_units
304		character(len=8) :: t_grid, t_grid_read
305		character(len=64) :: ax_name(5) ! The axis names of a variable
306
307	4	use_lath = .false. ; use_lonh = .false.
308	4	use_latq = .false. ; use_lonq = .false.
309	4	use_time = .false. ; use_periodic = .false.
310	4	use_layer = .false. ; use_int = .false.
311	4	num_extra_dims = 0
312	4	if (present(extra_axes)) then
313	2	num_extra_dims = size(extra_axes)
314	2	if (num_extra_dims > 0) then
315	12	allocate(use_extra_axis(num_extra_dims)) ; use_extra_axis = .false.
316	12	allocate(more_axes(num_extra_dims))
317		endif
318		endif
319
320	4	thread = SINGLE_FILE
321	4	if (PRESENT(threading)) thread = threading
322
323	4	domain_set = .false.
324	4	if (present(G)) then
325	3	domain_set = .true. ; Domain => G%Domain
326	3	gridLatT => G%gridLatT ; gridLatB => G%gridLatB
327	3	gridLonT => G%gridLonT ; gridLonB => G%gridLonB
328	3	x_axis_units = G%x_axis_units ; y_axis_units = G%y_axis_units
329	3	isg = G%isg ; ieg = G%ieg ; jsg = G%jsg ; jeg = G%jeg
330	3	IsgB = G%IsgB ; IegB = G%IegB ; JsgB = G%JsgB ; JegB = G%JegB
331	1	elseif (present(dG)) then
332	0	domain_set = .true. ; Domain => dG%Domain
333	0	gridLatT => dG%gridLatT ; gridLatB => dG%gridLatB
334	0	gridLonT => dG%gridLonT ; gridLonB => dG%gridLonB
335	0	x_axis_units = dG%x_axis_units ; y_axis_units = dG%y_axis_units
336	0	isg = dG%isg ; ieg = dG%ieg ; jsg = dG%jsg ; jeg = dG%jeg
337	0	IsgB = dG%IsgB ; IegB = dG%IegB ; JsgB = dG%JsgB ; JegB = dG%JegB
338		endif
339
340	4	one_file = .true.
341	4	if (domain_set) one_file = (thread == SINGLE_FILE)
342
343	4	if (one_file) then
344	2	if (domain_set) then
345		call IO_handle%open(filename, action=OVERWRITE_FILE, &
346	1	MOM_domain=domain, threading=thread, fileset=SINGLE_FILE)
347		else
348		call IO_handle%open(filename, action=OVERWRITE_FILE, threading=thread, &
349	1	fileset=SINGLE_FILE)
350		endif
351		else
352		call IO_handle%open(filename, action=OVERWRITE_FILE, MOM_domain=Domain, &
353	2	threading=thread, fileset=thread)
354		endif
355
356		! Define the coordinates.
357	87	do k=1,novars
358	83	position = vars(k)%position
359	83	if (position == -1) position = position_from_horgrid(vars(k)%hor_grid)
360	37	select case (position)
361	37	case (CENTER) ; use_lath = .true. ; use_lonh = .true.
362	2	case (CORNER) ; use_latq = .true. ; use_lonq = .true.
363	10	case (EAST_FACE) ; use_lath = .true. ; use_lonq = .true.
364	34	case (NORTH_FACE) ; use_latq = .true. ; use_lonh = .true.
365		case (0) ! Do nothing.
366		case default
367	83	call MOM_error(WARNING, "MOM_io create_file: "//trim(vars(k)%name)//" has an unrecognized value of postion")
368		end select
369	110	select case (vars(k)%z_grid)
370	27	case ('L') ; use_layer = .true.
371	56	case ('i') ; use_int = .true.
372		case ('1') ! Do nothing.
373		case default
374		call MOM_error(FATAL, "MOM_io create_file: "//trim(vars(k)%name)//&
375	83	" has unrecognized z_grid "//trim(vars(k)%z_grid))
376		end select
377	83	t_grid = adjustl(vars(k)%t_grid)
378	81	select case (t_grid(1:1))
379	81	case ('s', 'a', 'm') ; use_time = .true.
380	0	case ('p') ; use_periodic = .true.
381	0	if (len_trim(t_grid(2:8)) <= 0) call MOM_error(FATAL, &
382		"MOM_io create_file: No periodic axis length was specified in "//&
383		trim(vars(k)%t_grid) // " in the periodic axes of variable "//&
384	0	trim(vars(k)%name)//" in file "//trim(filename))
385	0	var_periods = -9999999
386	0	t_grid_read = adjustl(t_grid(2:8))
387	0	read(t_grid_read,*) var_periods
388	0	if (var_periods == -9999999) call MOM_error(FATAL, &
389		"MOM_io create_file: Failed to read the number of periods from "//&
390		trim(vars(k)%t_grid) // " in the periodic axes of variable "//&
391	0	trim(vars(k)%name)//" in file "//trim(filename))
392	0	if (var_periods < 1) call MOM_error(FATAL, "MOM_io create_file: "//&
393		"variable "//trim(vars(k)%name)//" in file "//trim(filename)//&
394		" uses a periodic time axis, and must have a positive "//&
395	0	"value for the number of periods in "//vars(k)%t_grid )
396	0	if ((num_periods > 0) .and. (var_periods /= num_periods)) &
397		call MOM_error(FATAL, "MOM_io create_file: "//&
398		"Only one value of the number of periods can be used in the "//&
399		"create_file call for file "//trim(filename)//". The second is "//&
400	0	"variable "//trim(vars(k)%name)//" with t_grid "//vars(k)%t_grid )
401
402	2	num_periods = var_periods
403		case ('1') ! Do nothing.
404		case default
405		call MOM_error(WARNING, "MOM_io create_file: "//trim(vars(k)%name)//&
406	83	" has unrecognized t_grid "//trim(vars(k)%t_grid))
407		end select
408
409	502	do n=1,5 ; if (len_trim(vars(k)%dim_names(n)) > 0) then
410	0	dim_found = .false.
411	0	do m=1,num_extra_dims
412	0	if (lowercase(trim(vars(k)%dim_names(n))) == lowercase(trim(extra_axes(m)%name))) then
413	0	use_extra_axis(m) = .true.
414	0	dim_found = .true.
415	0	exit
416		endif
417		enddo
418	0	if (.not.dim_found) call MOM_error(FATAL, "Unable to find a match for dimension "//&
419	0	trim(vars(k)%dim_names(n))//" for variable "//trim(vars(k)%name)//" in file "//trim(filename))
420		endif ; enddo
421		enddo
422
423	4	if ((use_lath .or. use_lonh .or. use_latq .or. use_lonq)) then
424	2	if (.not.domain_set) call MOM_error(FATAL, "create_file: "//&
425	0	"An ocean_grid_type or dyn_horgrid_type is required to create a file with a horizontal coordinate.")
426
427	2	call get_domain_components(Domain, x_domain, y_domain)
428		endif
429	4	if ((use_layer .or. use_int) .and. .not.present(GV)) call MOM_error(FATAL, &
430	0	"create_file: A vertical grid type is required to create a file with a vertical coordinate.")
431
432	4	if (use_lath) &
433		axis_lath = IO_handle%register_axis("lath", units=y_axis_units, longname="Latitude", &
434	2	cartesian='Y', domain=y_domain, data=gridLatT(jsg:jeg))
435	4	if (use_lonh) &
436		axis_lonh = IO_handle%register_axis("lonh", units=x_axis_units, longname="Longitude", &
437	2	cartesian='X', domain=x_domain, data=gridLonT(isg:ieg))
438	4	if (use_latq) &
439		axis_latq = IO_handle%register_axis("latq", units=y_axis_units, longname="Latitude", &
440	2	cartesian='Y', domain=y_domain, data=gridLatB(JsgB:JegB), edge_axis=.true.)
441	4	if (use_lonq) &
442		axis_lonq = IO_handle%register_axis("lonq", units=x_axis_units, longname="Longitude", &
443	2	cartesian='X', domain=x_domain, data=gridLonB(IsgB:IegB), edge_axis=.true.)
444	4	if (use_layer) &
445		axis_layer = IO_handle%register_axis("Layer", units=trim(GV%zAxisUnits), &
446		longname="Layer "//trim(GV%zAxisLongName), cartesian='Z', &
447	4	sense=1, data=GV%sLayer(1:GV%ke))
448	4	if (use_int) &
449		axis_int = IO_handle%register_axis("Interface", units=trim(GV%zAxisUnits), &
450		longname="Interface "//trim(GV%zAxisLongName), cartesian='Z', &
451	4	sense=1, data=GV%sInterface(1:GV%ke+1))
452
453	4	if (use_time) then ; if (present(timeunit)) then
454		! Set appropriate units, depending on the value.
455	1	if (timeunit < 0.0) then
456	0	time_units = "days" ! The default value.
457	1	elseif ((timeunit >= 0.99) .and. (timeunit < 1.01)) then
458	0	time_units = "seconds"
459	1	elseif ((timeunit >= 3599.0) .and. (timeunit < 3601.0)) then
460	0	time_units = "hours"
461	1	elseif ((timeunit >= 86399.0) .and. (timeunit < 86401.0)) then
462	1	time_units = "days"
463	0	elseif ((timeunit >= 3.0e7) .and. (timeunit < 3.2e7)) then
464	0	time_units = "years"
465		else
466	0	write(time_units,'(es8.2," s")') timeunit
467		endif
468
469	1	axis_time = IO_handle%register_axis("Time", units=time_units, longname="Time", cartesian='T')
470		else
471	2	axis_time = IO_handle%register_axis("Time", units="days", longname="Time", cartesian='T')
472		endif ; endif
473
474	4	if (use_periodic) then
475	0	if (num_periods <= 1) call MOM_error(FATAL, "MOM_io create_file: "//&
476	0	"num_periods for file "//trim(filename)//" must be at least 1.")
477		! Define a periodic axis with unit labels.
478	0	allocate(axis_val(num_periods))
479	0	do k=1,num_periods ; axis_val(k) = real(k) ; enddo
480		axis_periodic = IO_handle%register_axis("Period", units="nondimensional", &
481	0	longname="Periods for cyclical variables", cartesian='T', data=axis_val)
482	0	deallocate(axis_val)
483		endif
484
485	14	do m=1,num_extra_dims ; if (use_extra_axis(m)) then
486	0	if (allocated(extra_axes(m)%ax_data)) then
487		more_axes(m) = IO_handle%register_axis(extra_axes(m)%name, units=extra_axes(m)%units, &
488		longname=extra_axes(m)%longname, cartesian=extra_axes(m)%cartesian, &
489	0	sense=extra_axes(m)%sense, data=extra_axes(m)%ax_data)
490	0	elseif (trim(extra_axes(m)%cartesian) == "T") then
491		more_axes(m) = IO_handle%register_axis(extra_axes(m)%name, units=extra_axes(m)%units, &
492	0	longname=extra_axes(m)%longname, cartesian=extra_axes(m)%cartesian)
493		else
494		! FMS requires that non-time axes have variables that label their values, even if they are trivial.
495	0	allocate (axis_val(extra_axes(m)%ax_size))
496	0	do k=1,extra_axes(m)%ax_size ; axis_val(k) = real(k) ; enddo
497		more_axes(m) = IO_handle%register_axis(extra_axes(m)%name, units=extra_axes(m)%units, &
498		longname=extra_axes(m)%longname, cartesian=extra_axes(m)%cartesian, &
499	0	sense=extra_axes(m)%sense, data=axis_val)
500	0	deallocate(axis_val)
501		endif
502		endif ; enddo
503
504	87	do k=1,novars
505	83	numaxes = 0
506	83	position = vars(k)%position
507	83	if (position == -1) position = position_from_horgrid(vars(k)%hor_grid)
508	37	select case (position)
509		case (CENTER)
510	37	numaxes = 2 ; axes(1) = axis_lonh ; axes(2) = axis_lath ; ax_name(1) = "lonh" ; ax_name(2) = "lath"
511		case (CORNER)
512	2	numaxes = 2 ; axes(1) = axis_lonq ; axes(2) = axis_latq ; ax_name(1) = "lonq" ; ax_name(2) = "latq"
513		case (EAST_FACE)
514	10	numaxes = 2 ; axes(1) = axis_lonq ; axes(2) = axis_lath ; ax_name(1) = "lonq" ; ax_name(2) = "lath"
515		case (NORTH_FACE)
516	34	numaxes = 2 ; axes(1) = axis_lonh ; axes(2) = axis_latq ; ax_name(1) = "lonh" ; ax_name(2) = "latq"
517		case (0) ! Do nothing.
518		case default
519		call MOM_error(WARNING, "MOM_io create_file: "//trim(vars(k)%name)//&
520	83	" has unrecognized position, hor_grid = "//trim(vars(k)%hor_grid))
521		end select
522	110	select case (vars(k)%z_grid)
523	27	case ('L') ; numaxes = numaxes+1 ; axes(numaxes) = axis_layer ; ax_name(numaxes) = "Layer"
524	56	case ('i') ; numaxes = numaxes+1 ; axes(numaxes) = axis_int ; ax_name(numaxes) = "Interface"
525		case ('1') ! Do nothing.
526		case default
527		call MOM_error(FATAL, "MOM_io create_file: "//trim(vars(k)%name)//&
528	83	" has unrecognized z_grid "//trim(vars(k)%z_grid))
529		end select
530
531	238	do n=1,numaxes
532	155	if ( (len_trim(vars(k)%dim_names(n)) > 0) .and. (trim(ax_name(n)) /= trim(vars(k)%dim_names(n))) ) &
533		call MOM_error(WARNING, "MOM_io create_file: dimension "//trim(ax_name(n))//&
534		" of variable "//trim(vars(k)%name)//" in "//trim(filename)//&
535	83	" is being set inconsistently as "//trim(vars(k)%dim_names(n)))
536		enddo
537	343	do n=numaxes+1,5 ; if (len_trim(vars(k)%dim_names(n)) > 0) then
538	0	dim_found = .false.
539	0	do m=1,num_extra_dims
540	0	if (lowercase(trim(vars(k)%dim_names(n))) == lowercase(trim(extra_axes(m)%name))) then
541	0	numaxes = numaxes+1 ; axes(numaxes) = more_axes(m)
542	0	exit
543		endif
544		enddo
545		endif ; enddo
546
547	83	t_grid = adjustl(vars(k)%t_grid)
548	81	select case (t_grid(1:1))
549	81	case ('s', 'a', 'm') ; numaxes = numaxes+1 ; axes(numaxes) = axis_time
550	2	case ('p') ; numaxes = numaxes+1 ; axes(numaxes) = axis_periodic
551		case ('1') ! Do nothing.
552		case default
553		call MOM_error(WARNING, "MOM_io create_file: "//trim(vars(k)%name)//&
554	83	" has unrecognized t_grid "//trim(vars(k)%t_grid))
555		end select
556
557	83	pack = 1
558	87	if (present(checksums)) then
559		fields(k) = IO_handle%register_field(axes(1:numaxes), vars(k)%name, vars(k)%units, &
560	63	vars(k)%longname, pack=pack, checksum=checksums(k,:), conversion=vars(k)%conversion)
561		else
562		fields(k) = IO_handle%register_field(axes(1:numaxes), vars(k)%name, vars(k)%units, &
563	20	vars(k)%longname, pack=pack, conversion=vars(k)%conversion)
564		endif
565		enddo
566
567	4	if (present(global_atts)) then
568	0	do n=1,size(global_atts)
569	0	if (allocated(global_atts(n)%name) .and. allocated(global_atts(n)%att_val)) &
570	0	call IO_handle%write_attribute(global_atts(n)%name, global_atts(n)%att_val)
571		enddo
572		endif
573
574		! Now write the variables with the axis label values
575	4	if (use_lath) call IO_handle%write_field(axis_lath)
576	4	if (use_latq) call IO_handle%write_field(axis_latq)
577	4	if (use_lonh) call IO_handle%write_field(axis_lonh)
578	4	if (use_lonq) call IO_handle%write_field(axis_lonq)
579	4	if (use_layer) call IO_handle%write_field(axis_layer)
580	4	if (use_int) call IO_handle%write_field(axis_int)
581	4	if (use_periodic) call IO_handle%write_field(axis_periodic)
582	14	do m=1,num_extra_dims ; if (use_extra_axis(m)) then
583	0	call IO_handle%write_field(more_axes(m))
584		endif ; enddo
585
586	4	if (num_extra_dims > 0) then
587	14	deallocate(use_extra_axis, more_axes)
588		endif
589	12	end subroutine create_MOM_file
590
591
592		!> `reopen_MOM_file` wrapper for the legacy file handle, `file_type`.
593		!! NOTE: This function may be removed in a future release.
594	0	subroutine reopen_file(IO_handle, filename, vars, novars, fields, threading, &
595	0	timeunit, G, dG, GV, extra_axes, global_atts)
596		type(file_type), intent(inout) :: IO_handle
597		!< Handle for a file or fileset that is to be opened or reopened for
598		!! writing
599		character(len=*), intent(in) :: filename
600		!< full path to the file to create
601		type(vardesc), intent(in) :: vars(:)
602		!< structures describing fields written to filename
603		integer, intent(in) :: novars
604		!< number of fields written to filename
605		type(fieldtype), intent(inout) :: fields(:)
606		!< array of fieldtypes for each variable
607		integer, optional, intent(in) :: threading
608		!< SINGLE_FILE or MULTIPLE
609		real, optional, intent(in) :: timeunit
610		!< length of the units for time [s]. The default value is 86400.0, for 1
611		!! day.
612		type(ocean_grid_type), optional, intent(in) :: G
613		!< ocean horizontal grid structure; G or dG is required if a new file uses
614		!! any horizontal grid axes.
615		type(dyn_horgrid_type), optional, intent(in) :: dG
616		!< dynamic horizontal grid structure; G or dG is required if a new file
617		!! uses any horizontal grid axes.
618		type(verticalGrid_type), optional, intent(in) :: GV
619		!< ocean vertical grid structure, which is required if a new file uses any
620		!! vertical grid axes.
621		type(axis_info), optional, intent(in) :: extra_axes(:)
622		!< Types with information about some axes that might be used in this file
623		type(attribute_info), optional, intent(in) :: global_atts(:)
624		!< Global attributes to write to this file
625
626	0	type(MOM_infra_file) :: mfile
627		!< Wrapper to MOM file
628	0	type(MOM_field), allocatable :: mfields(:)
629		!< Wrapper to MOM fields
630		integer :: i
631
632	0	mfile%handle_infra = IO_handle
633	0	allocate(mfields(size(fields)))
634
635		call reopen_MOM_file(mfile, filename, vars, novars, mfields, &
636		threading=threading, timeunit=timeunit, G=G, dG=dG, GV=GV, &
637	0	extra_axes=extra_axes, global_atts=global_atts)
638
639	0	IO_handle = mfile%handle_infra
640	0	call get_file_fields(IO_handle, fields)
641	0	end subroutine reopen_file
642
643
644		!> This routine opens an existing NetCDF file for output. If it
645		!! does not find the file, a new file is created. It also sets up
646		!! structures that describe this file and the variables that will
647		!! later be written to this file.
648	0	subroutine reopen_MOM_file(IO_handle, filename, vars, novars, fields, &
649	0	threading, timeunit, G, dG, GV, extra_axes, global_atts)
650		class(MOM_file), intent(inout) :: IO_handle !< Handle for a file or fileset that is to be
651		!! opened or reopened for writing
652		character(len=*), intent(in) :: filename !< full path to the file to create
653		type(vardesc), intent(in) :: vars(:) !< structures describing fields written to filename
654		integer, intent(in) :: novars !< number of fields written to filename
655		type(MOM_field), intent(inout) :: fields(:) !< array of fieldtypes for each variable
656		integer, optional, intent(in) :: threading !< SINGLE_FILE or MULTIPLE
657		real, optional, intent(in) :: timeunit !< length of the units for time [s]. The
658		!! default value is 86400.0, for 1 day.
659		type(ocean_grid_type), optional, intent(in) :: G !< ocean horizontal grid structure; G or dG
660		!! is required if a new file uses any
661		!! horizontal grid axes.
662		type(dyn_horgrid_type), optional, intent(in) :: dG !< dynamic horizontal grid structure; G or dG
663		!! is required if a new file uses any
664		!! horizontal grid axes.
665		type(verticalGrid_type), optional, intent(in) :: GV !< ocean vertical grid structure, which is
666		!! required if a new file uses any
667		!! vertical grid axes.
668		type(axis_info), optional, intent(in) :: extra_axes(:) !< Types with information about
669		!! some axes that might be used in this file
670		type(attribute_info), optional, intent(in) :: global_atts(:) !< Global attributes to
671		!! write to this file
672
673		type(MOM_domain_type), pointer :: Domain => NULL()
674		character(len=200) :: check_name, mesg
675		integer :: length, nvar, thread
676		logical :: exists, one_file, domain_set
677
678	0	thread = SINGLE_FILE
679	0	if (PRESENT(threading)) thread = threading
680
681		! For single-file IO, only the root PE is required to set up the fields.
682		! This permits calls by either the root PE or all PEs
683	0	if (.not. is_root_PE() .and. thread == SINGLE_FILE) return
684
685		! For multiple IO domains, we would need additional functionality:
686		! * Identify ranks as IO PEs
687		! * Determine the filename of
688		! Neither of these tasks should be handed by MOM6, so we cannot safely use
689		! this function. A framework-specific `inquire()` function is needed.
690		! Until it exists, we will disable this function.
691	0	if (thread == MULTIPLE) &
692		call MOM_error(FATAL, 'reopen_MOM_file does not yet support files with ' &
693	0	// 'multiple I/O domains.')
694
695	0	check_name = filename
696	0	length = len(trim(check_name))
697	0	if (check_name(length-2:length) /= ".nc") check_name = trim(check_name)//".nc"
698	0	if (thread /= SINGLE_FILE) check_name = trim(check_name)//".0000"
699
700	0	inquire(file=check_name,EXIST=exists)
701
702	0	if (.not.exists) then
703		call create_MOM_file(IO_handle, filename, vars, novars, fields, &
704		threading, timeunit, G=G, dG=dG, GV=GV, extra_axes=extra_axes, &
705	0	global_atts=global_atts)
706		else
707
708	0	domain_set = .false.
709	0	if (present(G)) then
710	0	domain_set = .true. ; Domain => G%Domain
711	0	elseif (present(dG)) then
712	0	domain_set = .true. ; Domain => dG%Domain
713		endif
714
715	0	one_file = .true.
716	0	if (domain_set) one_file = (thread == SINGLE_FILE)
717
718	0	if (one_file) then
719	0	call IO_handle%open(filename, APPEND_FILE, threading=thread)
720		else
721	0	call IO_handle%open(filename, APPEND_FILE, MOM_domain=Domain)
722		endif
723	0	if (.not. IO_handle%file_is_open()) return
724
725	0	call IO_handle%get_file_info(nvar=nvar)
726
727	0	if (nvar == -1) then
728	0	write (mesg,*) "Reopening file ",trim(filename)," apparently had ",nvar,&
729	0	" variables. Clobbering and creating file with ",novars," instead."
730	0	call MOM_error(WARNING,"MOM_io: "//mesg)
731		call create_MOM_file(IO_handle, filename, vars, novars, fields, &
732		threading, timeunit, G=G, dG=dG, GV=GV, extra_axes=extra_axes, &
733	0	global_atts=global_atts)
734	0	elseif (nvar /= novars) then
735	0	write (mesg,*) "Reopening file ",trim(filename)," with ",novars,&
736	0	" variables instead of ",nvar,"."
737	0	call MOM_error(FATAL,"MOM_io: "//mesg)
738		endif
739
740	0	if (nvar > 0) call IO_handle%get_file_fields(fields(1:nvar))
741		endif
742	0	end subroutine reopen_MOM_file
743
744
745		!> Return the index of sdtout if called from the root PE, or 0 for other PEs.
746	0	integer function stdout_if_root()
747	0	stdout_if_root = 0
748	0	if (is_root_PE()) stdout_if_root = stdout
749	0	end function stdout_if_root
750
751		!> This function determines how many time levels a variable has in a file.
752	0	function num_timelevels(filename, varname, min_dims) result(n_time)
753		character(len=*), intent(in) :: filename !< name of the file to read
754		character(len=*), intent(in) :: varname !< variable whose number of time levels
755		!! are to be returned
756		integer, optional, intent(in) :: min_dims !< The minimum number of dimensions a variable must have
757		!! if it has a time dimension. If the variable has 1 less
758		!! dimension than this, then 0 is returned.
759		integer :: n_time !< number of time levels varname has in filename
760
761		character(len=256) :: msg
762		integer :: ndims
763		integer :: sizes(8)
764
765	0	n_time = -1
766
767		! To do almost the same via MOM_io_infra calls, we could do the following:
768		! found = field_exists(filename, varname)
769		! if (found) then
770		! call open_file(ncid, filename, action=READONLY_FILE, form=NETCDF_FILE, threading=MULTIPLE)
771		! call get_file_info(ncid, ntime=n_time)
772		! endif
773		! However, this does not handle the case where the time axis for the variable is not the record
774		! axis and min_dims is not used.
775
776	0	call get_var_sizes(filename, varname, ndims, sizes, match_case=.false., caller="num_timelevels")
777
778	0	if (ndims > 0) n_time = sizes(ndims)
779
780	0	if (present(min_dims)) then
781	0	if (ndims < min_dims-1) then
782	0	write(msg, '(I0)') min_dims
783		call MOM_error(WARNING, "num_timelevels: variable "//trim(varname)//" in file "//&
784	0	trim(filename)//" has fewer than min_dims = "//trim(msg)//" dimensions.")
785	0	n_time = -1
786	0	elseif (ndims == min_dims - 1) then
787	0	n_time = 0
788		endif
789		endif
790
791	0	end function num_timelevels
792
793
794		!> get_var_sizes returns the number and size of dimensions associate with a variable in a file.
795		!! Usually only the root PE does the read, and then the information is broadcast
796	0	subroutine get_var_sizes(filename, varname, ndims, sizes, match_case, caller, all_read, dim_names, ncid_in)
797		character(len=*), intent(in) :: filename !< Name of the file to read, used here in messages
798		character(len=*), intent(in) :: varname !< The variable name, used here for messages
799		integer, intent(out) :: ndims !< The number of dimensions to the variable
800		integer, dimension(:), intent(out) :: sizes !< The dimension sizes, or 0 for extra values
801		logical, optional, intent(in) :: match_case !< If false, allow for variables name matches to be
802		!! case insensitive, but take a perfect match if
803		!! found. The default is true.
804		character(len=*), optional, intent(in) :: caller !< The name of a calling routine for use in error messages
805		logical, optional, intent(in) :: all_read !< If present and true, all PEs that call this
806		!! routine actually do the read, otherwise only
807		!! root PE reads and then it broadcasts the results.
808		character(len=*), dimension(:), &
809		optional, intent(out) :: dim_names !< The names of the dimensions for this variable
810		integer, optional, intent(in) :: ncid_in !< The netCDF ID of an open file. If absent, the
811		!! file is opened and closed within this routine.
812
813		logical :: do_read, do_broadcast
814	0	integer, allocatable :: size_msg(:) ! An array combining the number of dimensions and the sizes.
815		integer :: n, nval
816
817	0	do_read = is_root_pe()
818	0	if (present(all_read)) do_read = all_read .or. do_read
819	0	do_broadcast = .true. ; if (present(all_read)) do_broadcast = .not.all_read
820
821	0	if (do_read) call read_var_sizes(filename, varname, ndims, sizes, match_case, caller, dim_names, ncid_in)
822
823	0	if (do_broadcast) then
824		! Distribute the sizes from the root PE.
825	0	nval = size(sizes) + 1
826
827	0	allocate(size_msg(nval))
828	0	size_msg(1) = ndims
829	0	do n=2,nval ; size_msg(n) = sizes(n-1) ; enddo
830
831	0	call broadcast(size_msg, nval, blocking=.true.)
832
833	0	ndims = size_msg(1)
834	0	do n=2,nval ; sizes(n-1) = size_msg(n) ; enddo
835	0	deallocate(size_msg)
836
837	0	if (present(dim_names) .and. (ndims > 0)) then
838	0	nval = min(ndims, size(dim_names))
839	0	call broadcast(dim_names(1:nval), len(dim_names(1)), blocking=.true.)
840		endif
841		endif
842
843	0	end subroutine get_var_sizes
844
845		!> read_var_sizes returns the number and size of dimensions associated with a variable in a file.
846		!! If the variable is not in the file the returned sizes are all 0 and ndims is -1.
847		!! Every processor for which this is called does the reading.
848	0	subroutine read_var_sizes(filename, varname, ndims, sizes, match_case, caller, dim_names, ncid_in)
849		character(len=*), intent(in) :: filename !< Name of the file to read, used here in messages
850		character(len=*), intent(in) :: varname !< The variable name, used here for messages
851		integer, intent(out) :: ndims !< The number of dimensions to the variable
852		integer, dimension(:), intent(out) :: sizes !< The dimension sizes, or 0 for extra values
853		logical, optional, intent(in) :: match_case !< If false, allow for variables name matches to be
854		!! case insensitive, but take a perfect match if
855		!! found. The default is true.
856		character(len=*), &
857		optional, intent(in) :: caller !< The name of a calling routine for use in error messages
858		character(len=*), dimension(:), &
859		optional, intent(out) :: dim_names !< The names of the dimensions for this variable
860		integer, optional, intent(in) :: ncid_in !< The netCDF ID of an open file. If absent, the
861		!! file is opened and closed within this routine.
862
863		character(len=256) :: hdr, dimname
864	0	integer, allocatable :: dimids(:)
865		integer :: varid, ncid, n, status
866		logical :: success, found
867	0	hdr = "get_var_size: " ; if (present(caller)) hdr = trim(hdr)//": "
868	0	sizes(:) = 0 ; ndims = -1
869
870	0	if (present(ncid_in)) then
871	0	ncid = ncid_in
872		else
873	0	call open_file_to_read(filename, ncid, success=success)
874	0	if (.not.success) then
875	0	call MOM_error(WARNING, "Unsuccessfully attempted to open file "//trim(filename))
876	0	return
877		endif
878		endif
879
880		! Get the dimension sizes of the variable varname.
881	0	call get_varid(varname, ncid, filename, varid, match_case=match_case, found=found)
882	0	if (.not.found) then
883	0	call MOM_error(WARNING, "Could not find variable "//trim(varname)//" in file "//trim(filename))
884	0	return
885		endif
886
887	0	status = NF90_inquire_variable(ncid, varid, ndims=ndims)
888	0	if (status /= NF90_NOERR) then
889		call MOM_error(WARNING, trim(hdr) // trim(NF90_STRERROR(status)) //&
890	0	" Getting number of dimensions of "//trim(varname)//" in "//trim(filename))
891	0	return
892		endif
893	0	if (ndims < 1) return
894
895	0	allocate(dimids(ndims))
896	0	status = NF90_inquire_variable(ncid, varid, dimids=dimids(1:ndims))
897	0	if (status /= NF90_NOERR) then
898		call MOM_error(WARNING, trim(hdr) // trim(NF90_STRERROR(status)) //&
899	0	" Getting dimension IDs for "//trim(varname)//" in "//trim(filename))
900	0	deallocate(dimids) ; return
901		endif
902
903	0	do n = 1, min(ndims,size(sizes))
904	0	status = NF90_Inquire_Dimension(ncid, dimids(n), name=dimname, len=sizes(n))
905	0	if (status /= NF90_NOERR) call MOM_error(WARNING, trim(hdr) // trim(NF90_STRERROR(status)) //&
906	0	" Getting dimension length for "//trim(varname)//" in "//trim(filename))
907	0	if (present(dim_names)) then
908	0	if (n <= size(dim_names)) dim_names(n) = trim(dimname)
909		endif
910		enddo
911	0	deallocate(dimids)
912
913	0	if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)
914
915	0	end subroutine read_var_sizes
916
917		!> Read a real scalar variable from a netCDF file with the root PE, and broadcast the
918		!! results to all the other PEs.
919	0	subroutine read_variable_0d(filename, varname, var, ncid_in, scale)
920		character(len=*), intent(in) :: filename !< The name of the file to read
921		character(len=*), intent(in) :: varname !< The variable name of the data in the file
922		real, intent(inout) :: var !< The scalar into which to read the data in arbitrary units [A ~> a]
923		integer, optional, intent(in) :: ncid_in !< The netCDF ID of an open file. If absent, the
924		!! file is opened and closed within this routine
925		real, optional, intent(in) :: scale !< A scaling factor that the variable is multiplied by
926		!! before it is returned to convert from the units in the file
927		!! to the internal units for this variable [A a-1 ~> 1]
928
929		integer :: varid, ncid, rc
930		character(len=256) :: hdr
931	0	hdr = "read_variable_0d"
932
933	0	if (is_root_pe()) then
934	0	if (present(ncid_in)) then
935	0	ncid = ncid_in
936		else
937	0	call open_file_to_read(filename, ncid)
938		endif
939
940	0	call get_varid(varname, ncid, filename, varid, match_case=.false.)
941	0	if (varid < 0) call MOM_error(FATAL, "Unable to get netCDF varid for "//trim(varname)//&
942	0	" in "//trim(filename))
943	0	rc = NF90_get_var(ncid, varid, var)
944	0	if (rc /= NF90_NOERR) call MOM_error(FATAL, trim(hdr) // trim(NF90_STRERROR(rc)) //&
945	0	" Difficulties reading "//trim(varname)//" from "//trim(filename))
946
947	0	if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)
948
949	0	if (present(scale)) var = scale * var
950		endif
951
952	0	call broadcast(var, blocking=.true.)
953	0	end subroutine read_variable_0d
954
955		!> Read a 1-d real variable from a netCDF file with the root PE, and broadcast the
956		!! results to all the other PEs.
957	0	subroutine read_variable_1d(filename, varname, var, ncid_in, scale)
958		character(len=*), intent(in) :: filename !< The name of the file to read
959		character(len=*), intent(in) :: varname !< The variable name of the data in the file
960		real, dimension(:), intent(inout) :: var !< The 1-d array into which to read the data in arbitrary units [A ~> a]
961		integer, optional, intent(in) :: ncid_in !< The netCDF ID of an open file. If absent, the
962		!! file is opened and closed within this routine
963		real, optional, intent(in) :: scale !< A scaling factor that the variable is multiplied by
964		!! before it is returned to convert from the units in the file
965		!! to the internal units for this variable [A a-1 ~> 1]
966
967		integer :: varid, ncid, rc
968		character(len=256) :: hdr
969	0	hdr = "read_variable_1d"
970
971	0	if (is_root_pe()) then
972	0	if (present(ncid_in)) then
973	0	ncid = ncid_in
974		else
975	0	call open_file_to_read(filename, ncid)
976		endif
977
978	0	call get_varid(varname, ncid, filename, varid, match_case=.false.)
979	0	if (varid < 0) call MOM_error(FATAL, "Unable to get netCDF varid for "//trim(varname)//&
980	0	" in "//trim(filename))
981	0	rc = NF90_get_var(ncid, varid, var)
982	0	if (rc /= NF90_NOERR) call MOM_error(FATAL, trim(hdr) // trim(NF90_STRERROR(rc)) //&
983	0	" Difficulties reading "//trim(varname)//" from "//trim(filename))
984
985	0	if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)
986
987	0	if (present(scale)) then ; if (scale /= 1.0) then
988	0	var(:) = scale * var(:)
989		endif ; endif
990		endif
991
992	0	call broadcast(var, size(var), blocking=.true.)
993	0	end subroutine read_variable_1d
994
995		!> Read a integer scalar variable from a netCDF file with the root PE, and broadcast the
996		!! results to all the other PEs.
997	0	subroutine read_variable_0d_int(filename, varname, var, ncid_in)
998		character(len=*), intent(in) :: filename !< The name of the file to read
999		character(len=*), intent(in) :: varname !< The variable name of the data in the file
1000		integer, intent(inout) :: var !< The scalar into which to read the data
1001		integer, optional, intent(in) :: ncid_in !< The netCDF ID of an open file. If absent, the
1002		!! file is opened and closed within this routine.
1003
1004		integer :: varid, ncid, rc
1005		character(len=256) :: hdr
1006	0	hdr = "read_variable_0d_int"
1007
1008	0	if (is_root_pe()) then
1009	0	if (present(ncid_in)) then
1010	0	ncid = ncid_in
1011		else
1012	0	call open_file_to_read(filename, ncid)
1013		endif
1014
1015	0	call get_varid(varname, ncid, filename, varid, match_case=.false.)
1016	0	if (varid < 0) call MOM_error(FATAL, "Unable to get netCDF varid for "//trim(varname)//&
1017	0	" in "//trim(filename))
1018	0	rc = NF90_get_var(ncid, varid, var)
1019	0	if (rc /= NF90_NOERR) call MOM_error(FATAL, trim(hdr) // trim(NF90_STRERROR(rc)) //&
1020	0	" Difficulties reading "//trim(varname)//" from "//trim(filename))
1021
1022	0	if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)
1023		endif
1024
1025	0	call broadcast(var, blocking=.true.)
1026	0	end subroutine read_variable_0d_int
1027
1028		!> Read a 1-d integer variable from a netCDF file with the root PE, and broadcast the
1029		!! results to all the other PEs.
1030	0	subroutine read_variable_1d_int(filename, varname, var, ncid_in)
1031		character(len=*), intent(in) :: filename !< The name of the file to read
1032		character(len=*), intent(in) :: varname !< The variable name of the data in the file
1033		integer, dimension(:), intent(inout) :: var !< The 1-d array into which to read the data
1034		integer, optional, intent(in) :: ncid_in !< The netCDF ID of an open file. If absent, the
1035		!! file is opened and closed within this routine.
1036
1037		integer :: varid, ncid, rc
1038		character(len=256) :: hdr
1039	0	hdr = "read_variable_1d_int"
1040
1041	0	if (is_root_pe()) then
1042	0	if (present(ncid_in)) then
1043	0	ncid = ncid_in
1044		else
1045	0	call open_file_to_read(filename, ncid)
1046		endif
1047
1048	0	call get_varid(varname, ncid, filename, varid, match_case=.false.)
1049	0	if (varid < 0) call MOM_error(FATAL, "Unable to get netCDF varid for "//trim(varname)//&
1050	0	" in "//trim(filename))
1051	0	rc = NF90_get_var(ncid, varid, var)
1052	0	if (rc /= NF90_NOERR) call MOM_error(FATAL, trim(hdr) // trim(NF90_STRERROR(rc)) //&
1053	0	" Difficulties reading "//trim(varname)//" from "//trim(filename))
1054
1055	0	if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)
1056		endif
1057
1058	0	call broadcast(var, size(var), blocking=.true.)
1059	0	end subroutine read_variable_1d_int
1060
1061		!> Read a 2d array from a netCDF input file and save to a variable.
1062		!!
1063		!! Start and nread lenths may exceed var rank. This allows for reading slices
1064		!! of larger arrays.
1065		!!
1066		!! Previous versions of the model required a time axis on IO fields. This
1067		!! constraint was dropped in later versions. As a result, versions both with
1068		!! and without a time axis now exist. In order to support all such versions,
1069		!! we use a reshaped version of start and nread in order to read the variable
1070		!! as it exists in the file.
1071		!!
1072		!! Certain constraints are still applied to start and nread in order to ensure
1073		!! that varname is a valid 2d array, or contains valid 2d slices.
1074		!!
1075		!! I/O occurs only on the root PE, and data is broadcast to other ranks.
1076		!! Due to potentially large memory communication and storage, this subroutine
1077		!! should only be used when domain-decomposition is unavaialable.
1078	0	subroutine read_variable_2d(filename, varname, var, start, nread, ncid_in)
1079		character(len=*), intent(in) :: filename !< Name of file to be read
1080		character(len=*), intent(in) :: varname !< Name of variable to be read
1081		real, intent(out) :: var(:,:) !< Output array of variable [arbitrary]
1082		integer, optional, intent(in) :: start(:) !< Starting index on each axis.
1083		integer, optional, intent(in) :: nread(:) !< Number of values to be read along each axis
1084		integer, optional, intent(in) :: ncid_in !< netCDF ID of an opened file.
1085		!! If absent, the file is opened and closed within this routine.
1086
1087		integer :: ncid, varid
1088		integer :: field_ndims, dim_len
1089	0	integer, allocatable :: field_dimids(:), field_shape(:)
1090	0	integer, allocatable :: field_start(:), field_nread(:)
1091		integer :: i, rc
1092		character(len=*), parameter :: hdr = "read_variable_2d: "
1093
1094		! Validate shape of start and nread
1095	0	if (present(start)) then
1096	0	if (size(start) < 2) &
1097		call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) &
1098	0	// " start must have at least two dimensions.")
1099		endif
1100
1101	0	if (present(nread)) then
1102	0	if (size(nread) < 2) &
1103		call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) &
1104	0	// " nread must have at least two dimensions.")
1105
1106	0	if (any(nread(3:) > 1)) &
1107		call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) &
1108	0	// " nread may only read a single level in higher dimensions.")
1109		endif
1110
1111		! Since start and nread may be reshaped, we cannot rely on netCDF to ensure
1112		! that their lengths are equivalent, and must do it here.
1113	0	if (present(start) .and. present(nread)) then
1114	0	if (size(start) /= size(nread)) &
1115		call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) &
1116	0	// " start and nread must have the same length.")
1117		endif
1118
1119		! Open and read `varname` from `filename`
1120	0	if (is_root_pe()) then
1121	0	if (present(ncid_in)) then
1122	0	ncid = ncid_in
1123		else
1124	0	call open_file_to_Read(filename, ncid)
1125		endif
1126
1127	0	call get_varid(varname, ncid, filename, varid, match_case=.false.)
1128	0	if (varid < 0) call MOM_error(FATAL, "Unable to get netCDF varid for "//trim(varname)//&
1129	0	" in "//trim(filename))
1130
1131		! Query for the dimensionality of the input field
1132	0	rc = nf90_inquire_variable(ncid, varid, ndims=field_ndims)
1133	0	if (rc /= NF90_NOERR) call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) //&
1134	0	": Difficulties reading "//trim(varname)//" from "//trim(filename))
1135
1136		! Confirm that field is at least 2d
1137	0	if (field_ndims < 2) &
1138		call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) // " " // &
1139	0	trim(varname) // " from " // trim(filename) // " is not a 2D field.")
1140
1141		! If start and nread are present, then reshape them to match field dims
1142	0	if (present(start) .or. present(nread)) then
1143	0	allocate(field_shape(field_ndims))
1144	0	allocate(field_dimids(field_ndims))
1145
1146	0	rc = nf90_inquire_variable(ncid, varid, dimids=field_dimids)
1147	0	if (rc /= NF90_NOERR) call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) //&
1148	0	": Difficulties reading "//trim(varname)//" from "//trim(filename))
1149
1150	0	do i = 1, field_ndims
1151	0	rc = nf90_inquire_dimension(ncid, field_dimids(i), len=dim_len)
1152	0	if (rc /= NF90_NOERR) &
1153		call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) &
1154	0	// ": Difficulties reading dimensions from " // trim(filename))
1155	0	field_shape(i) = dim_len
1156		enddo
1157
1158		! Reshape start(:) and nreads(:) in case ranks differ
1159	0	allocate(field_start(field_ndims))
1160	0	field_start(:) = 1
1161	0	if (present(start)) then
1162	0	dim_len = min(size(start), size(field_start))
1163	0	field_start(:dim_len) = start(:dim_len)
1164		endif
1165
1166	0	allocate(field_nread(field_ndims))
1167	0	field_nread(:2) = field_shape(:2)
1168	0	field_nread(3:) = 1
1169	0	if (present(nread)) field_nread(:2) = nread(:2)
1170
1171	0	rc = nf90_get_var(ncid, varid, var, field_start, field_nread)
1172
1173	0	deallocate(field_start)
1174	0	deallocate(field_nread)
1175	0	deallocate(field_shape)
1176	0	deallocate(field_dimids)
1177		else
1178	0	rc = nf90_get_var(ncid, varid, var)
1179		endif
1180
1181	0	if (rc /= NF90_NOERR) call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) //&
1182	0	" Difficulties reading "//trim(varname)//" from "//trim(filename))
1183
1184	0	if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)
1185		endif
1186
1187	0	call broadcast(var, size(var), blocking=.true.)
1188	0	end subroutine read_variable_2d
1189
1190
1191	0	subroutine read_variable_3d(filename, varname, var, start, nread, ncid_in)
1192		character(len=*), intent(in) :: filename !< Name of file to be read
1193		character(len=*), intent(in) :: varname !< Name of variable to be read
1194		real, intent(out) :: var(:,:,:) !< Output array of variable [arbitrary]
1195		integer, optional, intent(in) :: start(:) !< Starting index on each axis.
1196		integer, optional, intent(in) :: nread(:) !< Number of values to be read along each axis
1197		integer, optional, intent(in) :: ncid_in !< netCDF ID of an opened file.
1198		!! If absent, the file is opened and closed within this routine.
1199
1200		integer :: ncid, varid
1201		integer :: field_ndims, dim_len
1202	0	integer, allocatable :: field_dimids(:), field_shape(:)
1203	0	integer, allocatable :: field_start(:), field_nread(:)
1204		integer :: i, rc
1205		character(len=*), parameter :: hdr = "read_variable_3d: "
1206
1207		! Validate shape of start and nread
1208	0	if (present(start)) then
1209	0	if (size(start) < 2) &
1210		call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) &
1211	0	// " start must have at least two dimensions.")
1212		endif
1213
1214	0	if (present(nread)) then
1215	0	if (size(nread) < 2) &
1216		call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) &
1217	0	// " nread must have at least two dimensions.")
1218
1219	0	if (any(nread(3:) > 1)) &
1220		call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) &
1221	0	// " nread may only read a single level in higher dimensions.")
1222		endif
1223
1224		! Since start and nread may be reshaped, we cannot rely on netCDF to ensure
1225		! that their lengths are equivalent, and must do it here.
1226	0	if (present(start) .and. present(nread)) then
1227	0	if (size(start) /= size(nread)) &
1228		call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) &
1229	0	// " start and nread must have the same length.")
1230		endif
1231
1232		! Open and read `varname` from `filename`
1233	0	if (is_root_pe()) then
1234	0	if (present(ncid_in)) then
1235	0	ncid = ncid_in
1236		else
1237	0	call open_file_to_Read(filename, ncid)
1238		endif
1239
1240	0	call get_varid(varname, ncid, filename, varid, match_case=.false.)
1241	0	if (varid < 0) call MOM_error(FATAL, "Unable to get netCDF varid for "//trim(varname)//&
1242	0	" in "//trim(filename))
1243
1244		! Query for the dimensionality of the input field
1245	0	rc = nf90_inquire_variable(ncid, varid, ndims=field_ndims)
1246	0	if (rc /= NF90_NOERR) call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) //&
1247	0	": Difficulties reading "//trim(varname)//" from "//trim(filename))
1248
1249		! Confirm that field is at least 2d
1250	0	if (field_ndims < 2) &
1251		call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) // " " // &
1252	0	trim(varname) // " from " // trim(filename) // " is not a 2D field.")
1253
1254		! If start and nread are present, then reshape them to match field dims
1255	0	if (present(start) .or. present(nread)) then
1256	0	allocate(field_shape(field_ndims))
1257	0	allocate(field_dimids(field_ndims))
1258
1259	0	rc = nf90_inquire_variable(ncid, varid, dimids=field_dimids)
1260	0	if (rc /= NF90_NOERR) call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) //&
1261	0	": Difficulties reading "//trim(varname)//" from "//trim(filename))
1262
1263	0	do i = 1, field_ndims
1264	0	rc = nf90_inquire_dimension(ncid, field_dimids(i), len=dim_len)
1265	0	if (rc /= NF90_NOERR) &
1266		call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) &
1267	0	// ": Difficulties reading dimensions from " // trim(filename))
1268	0	field_shape(i) = dim_len
1269		enddo
1270
1271		! Reshape start(:) and nreads(:) in case ranks differ
1272	0	allocate(field_start(field_ndims))
1273	0	field_start(:) = 1
1274	0	if (present(start)) then
1275	0	dim_len = min(size(start), size(field_start))
1276	0	field_start(:dim_len) = start(:dim_len)
1277		endif
1278
1279	0	allocate(field_nread(field_ndims))
1280	0	field_nread(:3) = field_shape(:3)
1281		!field_nread(3:) = 1
1282	0	if (present(nread)) field_nread(:3) = nread(:3)
1283
1284	0	rc = nf90_get_var(ncid, varid, var, field_start, field_nread)
1285
1286	0	deallocate(field_start)
1287	0	deallocate(field_nread)
1288	0	deallocate(field_shape)
1289	0	deallocate(field_dimids)
1290		else
1291	0	rc = nf90_get_var(ncid, varid, var)
1292		endif
1293
1294	0	if (rc /= NF90_NOERR) call MOM_error(FATAL, hdr // trim(nf90_strerror(rc)) //&
1295	0	" Difficulties reading "//trim(varname)//" from "//trim(filename))
1296
1297	0	if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)
1298		endif
1299
1300	0	call broadcast(var, size(var), blocking=.true.)
1301	0	end subroutine read_variable_3d
1302
1303		!> Read a character-string global or variable attribute
1304	0	subroutine read_attribute_str(filename, attname, att_val, varname, found, all_read, ncid_in)
1305		character(len=*), intent(in) :: filename !< Name of the file to read
1306		character(len=*), intent(in) :: attname !< Name of the attribute to read
1307		character(:), allocatable, intent(out) :: att_val !< The value of the attribute
1308		character(len=*), optional, intent(in) :: varname !< The name of the variable whose attribute will
1309		!! be read. If missing, read a global attribute.
1310		logical, optional, intent(out) :: found !< Returns true if the attribute is found
1311		logical, optional, intent(in) :: all_read !< If present and true, all PEs that call this
1312		!! routine actually do the read, otherwise only
1313		!! root PE reads and then broadcasts the results.
1314		integer, optional, intent(in) :: ncid_in !< The netCDF ID of an open file. If absent, the
1315		!! file is opened and closed within this routine.
1316
1317		logical :: do_read, do_broadcast
1318		integer :: rc, ncid, varid, att_type, att_len, info(2)
1319		character(len=256) :: hdr, att_str
1320	0	character(len=:), dimension(:), allocatable :: tmp_str
1321	0	hdr = "read_attribute_str"
1322	0	att_len = 0
1323
1324	0	do_read = is_root_pe() ; if (present(all_read)) do_read = all_read .or. do_read
1325	0	do_broadcast = .true. ; if (present(all_read)) do_broadcast = .not.all_read
1326
1327	0	if (do_read) then
1328	0	if (present(ncid_in)) then
1329	0	ncid = ncid_in
1330		else
1331	0	call open_file_to_read(filename, ncid, success=found)
1332	0	if (present(found)) then ; if (.not.found) do_read = .false. ; endif
1333		endif
1334		endif
1335
1336	0	if (do_read) then
1337	0	rc = NF90_ENOTATT ; att_len = 0
1338	0	if (present(varname)) then ! Read a variable attribute
1339	0	call get_varid(varname, ncid, filename, varid, match_case=.false., found=found)
1340	0	att_str = "att "//trim(attname)//" for "//trim(varname)//" from "//trim(filename)
1341		else ! Read a global attribute
1342	0	varid = NF90_GLOBAL
1343	0	att_str = "global att "//trim(attname)//" from "//trim(filename)
1344		endif
1345	0	if ((varid > 0) .or. (varid == NF90_GLOBAL)) then ! The named variable does exist, and found would be true.
1346	0	rc = NF90_inquire_attribute(ncid, varid, attname, xtype=att_type, len=att_len)
1347	0	if ((.not. present(found)) .or. (rc /= NF90_ENOTATT)) then
1348	0	if ((rc /= NF90_NOERR) .and. (rc /= NF90_ENOTATT)) &
1349	0	call MOM_error(FATAL, trim(hdr) // trim(NF90_STRERROR(rc)) //" Error getting info for "//trim(att_str))
1350	0	if (att_type /= NF90_CHAR) &
1351	0	call MOM_error(FATAL, trim(hdr)//": Attribute data type is not a char for "//trim(att_str))
1352		! if (att_len > len(att_val)) &
1353		! call MOM_error(FATAL, trim(hdr)//": Insufficiently long string passed in to read "//trim(att_str))
1354	0	allocate(character(att_len) :: att_val)
1355
1356	0	if (rc == NF90_NOERR) then
1357	0	rc = NF90_get_att(ncid, varid, attname, att_val)
1358	0	if ((rc /= NF90_NOERR) .and. (rc /= NF90_ENOTATT)) &
1359	0	call MOM_error(FATAL, trim(hdr) // trim(NF90_STRERROR(rc)) //" Difficulties reading "//trim(att_str))
1360		endif
1361		endif
1362		endif
1363	0	if (present(found)) found = (rc == NF90_NOERR)
1364
1365	0	if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)
1366		endif
1367
1368	0	if (do_broadcast) then
1369		! Communicate the string length
1370	0	info(1) = att_len ; info(2) = 0 ; if (do_read .and. found) info(2) = 1
1371	0	call broadcast(info, 2, blocking=.true.)
1372	0	if (present(found)) then
1373	0	found = (info(2) /= 0)
1374	0	if (.not. found) return
1375		endif
1376	0	att_len = info(1)
1377
1378	0	if (att_len > 0) then
1379		! These extra copies are here because broadcast only supports arrays of strings.
1380	0	allocate(character(att_len) :: tmp_str(1))
1381	0	if (.not.do_read) allocate(character(att_len) :: att_val)
1382	0	if (do_read) tmp_str(1) = att_val
1383	0	call broadcast(tmp_str, att_len, blocking=.true.)
1384	0	att_val = tmp_str(1)
1385	0	elseif (.not.allocated(att_val)) then
1386	0	allocate(character(4) :: att_val) ; att_val = ''
1387		endif
1388	0	elseif (.not.allocated(att_val)) then
1389	0	allocate(character(4) :: att_val) ; att_val = ''
1390		endif
1391	0	end subroutine read_attribute_str
1392
1393
1394		!> Read a 32-bit integer global or variable attribute
1395	0	subroutine read_attribute_int32(filename, attname, att_val, varname, found, all_read, ncid_in)
1396		character(len=*), intent(in) :: filename !< Name of the file to read
1397		character(len=*), intent(in) :: attname !< Name of the attribute to read
1398		integer(kind=int32), intent(out) :: att_val !< The value of the attribute
1399		character(len=*), optional, intent(in) :: varname !< The name of the variable whose attribute will
1400		!! be read. If missing, read a global attribute.
1401		logical, optional, intent(out) :: found !< Returns true if the attribute is found
1402		logical, optional, intent(in) :: all_read !< If present and true, all PEs that call this
1403		!! routine actually do the read, otherwise only
1404		!! root PE reads and then broadcasts the results.
1405		integer, optional, intent(in) :: ncid_in !< The netCDF ID of an open file. If absent, the
1406		!! file is opened and closed within this routine.
1407
1408		logical :: do_read, do_broadcast
1409		integer :: rc, ncid, varid, is_found
1410		character(len=256) :: hdr
1411	0	hdr = "read_attribute_int32"
1412	0	att_val = 0
1413
1414	0	do_read = is_root_pe() ; if (present(all_read)) do_read = all_read .or. do_read
1415	0	do_broadcast = .true. ; if (present(all_read)) do_broadcast = .not.all_read
1416
1417	0	if (do_read) then
1418	0	if (present(ncid_in)) then
1419	0	ncid = ncid_in
1420		else
1421	0	call open_file_to_read(filename, ncid, success=found)
1422	0	if (present(found)) then ; if (.not.found) do_read = .false. ; endif
1423		endif
1424		endif
1425
1426	0	if (do_read) then
1427	0	rc = NF90_ENOTATT
1428	0	if (present(varname)) then ! Read a variable attribute
1429	0	call get_varid(varname, ncid, filename, varid, match_case=.false., found=found)
1430	0	if (varid >= 0) then ! The named variable does exist, and found would be true.
1431	0	rc = NF90_get_att(ncid, varid, attname, att_val)
1432	0	if ((rc /= NF90_NOERR) .and. (rc /= NF90_ENOTATT)) &
1433		call MOM_error(FATAL, trim(hdr) // trim(NF90_STRERROR(rc)) //" Difficulties reading att "//&
1434	0	trim(attname)//" for "//trim(varname)//" from "//trim(filename))
1435		endif
1436		else ! Read a global attribute
1437	0	rc = NF90_get_att(ncid, NF90_GLOBAL, attname, att_val)
1438	0	if ((rc /= NF90_NOERR) .and. (rc /= NF90_ENOTATT)) &
1439		call MOM_error(FATAL, trim(hdr) // trim(NF90_STRERROR(rc)) //&
1440	0	" Difficulties reading global att "//trim(attname)//" from "//trim(filename))
1441		endif
1442	0	if (present(found)) found = (rc == NF90_NOERR)
1443
1444	0	if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)
1445		endif
1446
1447	0	if (do_broadcast) then
1448	0	if (present(found)) then
1449	0	is_found = 0 ; if (is_root_pe() .and. found) is_found = 1
1450	0	call broadcast(is_found, blocking=.false.)
1451		endif
1452	0	call broadcast(att_val, blocking=.true.)
1453	0	if (present(found)) found = (is_found /= 0)
1454		endif
1455
1456	0	end subroutine read_attribute_int32
1457
1458
1459		!> Read a 64-bit integer global or variable attribute
1460	0	subroutine read_attribute_int64(filename, attname, att_val, varname, found, all_read, ncid_in)
1461		character(len=*), intent(in) :: filename !< Name of the file to read
1462		character(len=*), intent(in) :: attname !< Name of the attribute to read
1463		integer(kind=int64), intent(out) :: att_val !< The value of the attribute
1464		character(len=*), optional, intent(in) :: varname !< The name of the variable whose attribute will
1465		!! be read. If missing, read a global attribute.
1466		logical, optional, intent(out) :: found !< Returns true if the attribute is found
1467		logical, optional, intent(in) :: all_read !< If present and true, all PEs that call this
1468		!! routine actually do the read, otherwise only
1469		!! root PE reads and then broadcasts the results.
1470		integer, optional, intent(in) :: ncid_in !< The netCDF ID of an open file. If absent, the
1471		!! file is opened and closed within this routine.
1472
1473		logical :: do_read, do_broadcast
1474		integer :: rc, ncid, varid, is_found
1475		character(len=256) :: hdr
1476	0	hdr = "read_attribute_int64"
1477	0	att_val = 0
1478
1479	0	do_read = is_root_pe() ; if (present(all_read)) do_read = all_read .or. do_read
1480	0	do_broadcast = .true. ; if (present(all_read)) do_broadcast = .not.all_read
1481
1482	0	if (do_read) then
1483	0	if (present(ncid_in)) then
1484	0	ncid = ncid_in
1485		else
1486	0	call open_file_to_read(filename, ncid, success=found)
1487	0	if (present(found)) then ; if (.not.found) do_read = .false. ; endif
1488		endif
1489		endif
1490
1491	0	if (do_read) then
1492	0	rc = NF90_ENOTATT
1493	0	if (present(varname)) then ! Read a variable attribute
1494	0	call get_varid(varname, ncid, filename, varid, match_case=.false., found=found)
1495	0	if (varid >= 0) then ! The named variable does exist, and found would be true.
1496	0	rc = NF90_get_att(ncid, varid, attname, att_val)
1497	0	if ((rc /= NF90_NOERR) .and. (rc /= NF90_ENOTATT)) &
1498		call MOM_error(FATAL, trim(hdr) // trim(NF90_STRERROR(rc)) //" Difficulties reading att "//&
1499	0	trim(attname)//" for "//trim(varname)//" from "//trim(filename))
1500		endif
1501		else ! Read a global attribute
1502	0	rc = NF90_get_att(ncid, NF90_GLOBAL, attname, att_val)
1503	0	if ((rc /= NF90_NOERR) .and. (rc /= NF90_ENOTATT)) &
1504		call MOM_error(FATAL, trim(hdr) // trim(NF90_STRERROR(rc)) //&
1505	0	" Difficulties reading global att "//trim(attname)//" from "//trim(filename))
1506		endif
1507	0	if (present(found)) found = (rc == NF90_NOERR)
1508
1509	0	rc = NF90_close(ncid)
1510		endif
1511
1512	0	if (do_broadcast) then
1513	0	if (present(found)) then
1514	0	is_found = 0 ; if (is_root_pe() .and. found) is_found = 1
1515	0	call broadcast(is_found, blocking=.false.)
1516		endif
1517	0	call broadcast(att_val, blocking=.true.)
1518	0	if (present(found)) found = (is_found /= 0)
1519		endif
1520
1521	0	end subroutine read_attribute_int64
1522
1523		!> Read a real global or variable attribute
1524	0	subroutine read_attribute_real(filename, attname, att_val, varname, found, all_read, ncid_in)
1525		character(len=*), intent(in) :: filename !< Name of the file to read
1526		character(len=*), intent(in) :: attname !< Name of the attribute to read
1527		real, intent(out) :: att_val !< The value of the attribute [arbitrary]
1528		character(len=*), optional, intent(in) :: varname !< The name of the variable whose attribute will
1529		!! be read. If missing, read a global attribute.
1530		logical, optional, intent(out) :: found !< Returns true if the attribute is found
1531		logical, optional, intent(in) :: all_read !< If present and true, all PEs that call this
1532		!! routine actually do the read, otherwise only
1533		!! root PE reads and then broadcasts the results.
1534		integer, optional, intent(in) :: ncid_in !< The netCDF ID of an open file. If absent, the
1535		!! file is opened and closed within this routine.
1536
1537		logical :: do_read, do_broadcast
1538		integer :: rc, ncid, varid, is_found
1539		character(len=256) :: hdr
1540	0	hdr = "read_attribute_real"
1541	0	att_val = 0.0
1542
1543	0	do_read = is_root_pe() ; if (present(all_read)) do_read = all_read .or. do_read
1544	0	do_broadcast = .true. ; if (present(all_read)) do_broadcast = .not.all_read
1545
1546	0	if (do_read) then
1547	0	if (present(ncid_in)) then
1548	0	ncid = ncid_in
1549		else
1550	0	call open_file_to_read(filename, ncid, success=found)
1551	0	if (present(found)) then ; if (.not.found) do_read = .false. ; endif
1552		endif
1553		endif
1554
1555	0	if (do_read) then
1556	0	rc = NF90_ENOTATT
1557	0	if (present(varname)) then ! Read a variable attribute
1558	0	call get_varid(varname, ncid, filename, varid, match_case=.false., found=found)
1559	0	if (varid >= 0) then ! The named variable does exist, and found would be true.
1560	0	rc = NF90_get_att(ncid, varid, attname, att_val)
1561	0	if ((rc /= NF90_NOERR) .and. (rc /= NF90_ENOTATT)) &
1562		call MOM_error(FATAL, trim(hdr) // trim(NF90_STRERROR(rc)) //" Difficulties reading att "//&
1563	0	trim(attname)//" for "//trim(varname)//" from "//trim(filename))
1564		endif
1565		else ! Read a global attribute
1566	0	rc = NF90_get_att(ncid, NF90_GLOBAL, attname, att_val)
1567	0	if ((rc /= NF90_NOERR) .and. (rc /= NF90_ENOTATT)) &
1568		call MOM_error(FATAL, trim(hdr) // trim(NF90_STRERROR(rc)) //&
1569	0	" Difficulties reading global att "//trim(attname)//" from "//trim(filename))
1570		endif
1571	0	if (present(found)) found = (rc == NF90_NOERR)
1572
1573	0	if (.not.present(ncid_in)) call close_file_to_read(ncid, filename)
1574		endif
1575
1576	0	if (do_broadcast) then
1577	0	if (present(found)) then
1578	0	is_found = 0 ; if (is_root_pe() .and. found) is_found = 1
1579	0	call broadcast(is_found, blocking=.false.)
1580		endif
1581	0	call broadcast(att_val, blocking=.true.)
1582	0	if (present(found)) found = (is_found /= 0)
1583		endif
1584
1585	0	end subroutine read_attribute_real
1586
1587		!> Open a netcdf file for reading, with error handling
1588	0	subroutine open_file_to_read(filename, ncid, success)
1589		character(len=*), intent(in) :: filename !< path and name of the file to open for reading
1590		integer, intent(out) :: ncid !< The netcdf handle for the file
1591		logical, optional, intent(out) :: success !< Returns true if the file was opened, or if this
1592		!! argument is not present, failure is fatal error.
1593		! Local variables
1594		integer rc
1595
1596	0	rc = NF90_open(trim(filename), NF90_NOWRITE, ncid)
1597	0	if (present(success)) then
1598	0	success = (rc == NF90_NOERR)
1599	0	elseif (rc /= NF90_NOERR) then
1600	0	call MOM_error(FATAL, "Difficulties opening "//trim(filename)//" - "//trim(NF90_STRERROR(rc)) )
1601		endif
1602
1603	0	end subroutine open_file_to_read
1604
1605		!> Close a netcdf file that had been opened for reading, with error handling
1606	0	subroutine close_file_to_read(ncid, filename)
1607		integer, intent(inout) :: ncid !< The netcdf handle for the file to close
1608		character(len=*), optional, intent(in) :: filename !< path and name of the file to close
1609		integer :: rc
1610	0	if (ncid >= 0) then
1611	0	rc = NF90_close(ncid)
1612	0	if (present(filename) .and. (rc /= NF90_NOERR)) then
1613	0	call MOM_error(WARNING, "Difficulties closing "//trim(filename)//": "//trim(NF90_STRERROR(rc)))
1614	0	elseif (rc /= NF90_NOERR) then
1615	0	call MOM_error(WARNING, "Difficulties closing file: "//trim(NF90_STRERROR(rc)))
1616		endif
1617		endif
1618	0	ncid = -1
1619	0	end subroutine close_file_to_read
1620
1621		!> get_varid finds the netcdf handle for the potentially case-insensitive variable name in a file
1622	0	subroutine get_varid(varname, ncid, filename, varid, match_case, found)
1623		character(len=*), intent(in) :: varname !< The name of the variable that is being sought
1624		integer, intent(in) :: ncid !< The open netcdf handle for the file
1625		character(len=*), intent(in) :: filename !< name of the file to read, used here in messages
1626		integer, intent(out) :: varid !< The netcdf handle for the variable
1627		logical, optional, intent(in) :: match_case !< If false, allow for variables name matches to be
1628		!! case insensitive, but take a perfect match if
1629		!! found. The default is true.
1630		logical, optional, intent(out) :: found !< Returns true if the attribute is found
1631
1632		logical :: var_found, insensitive
1633		character(len=256) :: name
1634	0	integer, allocatable :: varids(:)
1635		integer :: nvars, status, n
1636
1637	0	varid = -1
1638	0	var_found = .false.
1639	0	insensitive = .false. ; if (present(match_case)) insensitive = .not.match_case
1640
1641	0	if (insensitive) then
1642		! This code ounddoes a case-insensitive search for a variable in the file.
1643	0	status = NF90_inquire(ncid, nVariables=nvars)
1644	0	if (present(found) .and. ((status /= NF90_NOERR) .or. (nvars < 1))) then
1645	0	found = .false. ; return
1646	0	elseif (status /= NF90_NOERR) then
1647		call MOM_error(FATAL, "get_varid: Difficulties getting the number of variables in file "//&
1648	0	trim(filename)//" - "//trim(NF90_STRERROR(status)))
1649	0	elseif (nvars < 1) then
1650	0	call MOM_error(FATAL, "get_varid: There appear not to be any variables in "//trim(filename))
1651		endif
1652
1653	0	allocate(varids(nvars))
1654
1655	0	status = nf90_inq_varids(ncid, nvars, varids)
1656	0	if (status /= NF90_NOERR) then
1657		call MOM_error(WARNING, "get_varid: Difficulties getting the variable IDs in file "//&
1658	0	trim(filename)//" - "//trim(NF90_STRERROR(status)))
1659	0	nvars = -1 ! Full error handling will occur after the do-loop.
1660		endif
1661
1662	0	do n = 1,nvars
1663	0	status = nf90_inquire_variable(ncid, varids(n), name=name)
1664	0	if (status /= NF90_NOERR) then
1665		call MOM_error(WARNING, "get_varid: Difficulties getting a variable name in file "//&
1666	0	trim(filename)//" - "//trim(NF90_STRERROR(status)))
1667		endif
1668
1669	0	if (trim(lowercase(name)) == trim(lowercase(varname))) then
1670	0	if (var_found) then
1671		call MOM_error(WARNING, "get_varid: Two variables match the case-insensitive name "//&
1672	0	trim(varname)//" in file "//trim(filename))
1673		! Replace the first variable if the second one is a case-sensitive match
1674	0	if (trim(name) == trim(varname)) varid = varids(n)
1675		else
1676	0	varid = varids(n) ; var_found = .true.
1677		endif
1678		endif
1679		enddo
1680	0	if (present(found)) found = var_found
1681	0	if ((.not.var_found) .and. .not.present(found)) call MOM_error(FATAL, &
1682	0	"get_varid: variable "//trim(varname)//" was not found in file "//trim(filename))
1683
1684	0	deallocate(varids)
1685		else
1686	0	status = NF90_INQ_VARID(ncid, trim(varname), varid)
1687	0	if (present(found)) found = (status == NF90_NOERR)
1688	0	if ((status /= NF90_NOERR) .and. .not.present(found)) then
1689		call MOM_error(FATAL, "get_varid: Difficulties getting a variable id for "//&
1690	0	trim(varname)//" in file "//trim(filename)//" - "//trim(NF90_STRERROR(status)))
1691		endif
1692		endif
1693
1694	0	end subroutine get_varid
1695
1696		!> Verify that a file contains a named variable with the expected units.
1697	0	subroutine verify_variable_units(filename, varname, expected_units, msg, ierr, alt_units)
1698		character(len=*), intent(in) :: filename !< File name
1699		character(len=*), intent(in) :: varname !< Variable name
1700		character(len=*), intent(in) :: expected_units !< Expected units of variable
1701		character(len=*), intent(inout) :: msg !< Message to use for errors
1702		logical, intent(out) :: ierr !< True if an error occurs
1703		character(len=*), optional, intent(in) :: alt_units !< Alterate acceptable units of variable
1704
1705		! Local variables
1706		character (len=200) :: units
1707		logical :: units_correct, success
1708		integer :: i, ncid, status, vid
1709
1710	0	if (.not.is_root_pe()) then ! Only the root PE should do the verification.
1711	0	ierr = .false. ; msg = '' ; return
1712		endif
1713
1714	0	ierr = .true.
1715	0	call open_file_to_read(filename, ncid, success)
1716	0	if (.not.success) then
1717	0	msg = 'File not found: '//trim(filename)
1718	0	return
1719		endif
1720
1721	0	status = NF90_INQ_VARID(ncid, trim(varname), vid)
1722	0	if (status /= NF90_NOERR) then
1723	0	msg = 'Var not found: '//trim(varname)
1724		else
1725	0	status = NF90_GET_ATT(ncid, vid, "units", units)
1726	0	if (status /= NF90_NOERR) then
1727	0	msg = 'Attribute not found: units'
1728		else
1729		! NF90_GET_ATT can return attributes with null characters, which TRIM will not truncate.
1730		! This loop replaces any null characters with a space so that the subsequent check
1731		! between the read units and the expected units will pass
1732	0	do i=1,LEN_TRIM(units)
1733	0	if (units(i:i) == CHAR(0)) units(i:i) = " "
1734		enddo
1735
1736	0	units_correct = (trim(units) == trim(expected_units))
1737	0	if (present(alt_units)) then
1738	0	units_correct = units_correct .or. (trim(units) == trim(alt_units))
1739		endif
1740	0	if (units_correct) then
1741	0	ierr = .false.
1742	0	msg = ''
1743		else
1744	0	msg = 'Units incorrect: '//trim(units)//' /= '//trim(expected_units)
1745		endif
1746		endif
1747		endif
1748
1749	0	status = NF90_close(ncid)
1750
1751	0	end subroutine verify_variable_units
1752
1753		!> Returns a vardesc type whose elements have been filled with the provided
1754		!! fields. The argument name is required, while the others are optional and
1755		!! have default values that are empty strings or are appropriate for a 3-d
1756		!! tracer field at the tracer cell centers.
1757	91	function var_desc(name, units, longname, hor_grid, z_grid, t_grid, cmor_field_name, &
1758	91	cmor_units, cmor_longname, conversion, caller, position, dim_names, &
1759	91	extra_axes, fixed) result(vd)
1760		character(len=*), intent(in) :: name !< variable name
1761		character(len=*), optional, intent(in) :: units !< variable units
1762		character(len=*), optional, intent(in) :: longname !< variable long name
1763		character(len=*), optional, intent(in) :: hor_grid !< A character string indicating the horizontal
1764		!! position of this variable
1765		character(len=*), optional, intent(in) :: z_grid !< variable vertical staggering
1766		character(len=*), optional, intent(in) :: t_grid !< time description: s, p, or 1
1767		character(len=*), optional, intent(in) :: cmor_field_name !< CMOR name
1768		character(len=*), optional, intent(in) :: cmor_units !< CMOR physical dimensions of variable
1769		character(len=*), optional, intent(in) :: cmor_longname !< CMOR long name
1770		real , optional, intent(in) :: conversion !< for unit conversions, such as needed to
1771		!! convert from intensive to extensive
1772		!! [various] or [a A-1 ~> 1]
1773		character(len=*), optional, intent(in) :: caller !< The calling routine for error messages
1774		integer, optional, intent(in) :: position !< A coded integer indicating the horizontal position
1775		!! of this variable if it has such dimensions.
1776		!! Valid values include CORNER, CENTER, EAST_FACE
1777		!! NORTH_FACE, and 0 for no horizontal dimensions.
1778		character(len=*), dimension(:), &
1779		optional, intent(in) :: dim_names !< The names of the dimensions of this variable
1780		type(axis_info), dimension(:), &
1781		optional, intent(in) :: extra_axes !< dimensions other than space-time
1782		logical, optional, intent(in) :: fixed !< If true, this does not evolve with time
1783		type(vardesc) :: vd !< vardesc type that is created
1784
1785		character(len=120) :: cllr
1786	91	cllr = "var_desc"
1787	91	if (present(caller)) cllr = trim(caller)
1788
1789	91	call safe_string_copy(name, vd%name, "vd%name", cllr)
1790
1791	91	vd%longname = "" ; vd%units = ""
1792	91	vd%hor_grid = 'h' ; vd%position = CENTER ; vd%z_grid = 'L' ; vd%t_grid = 's'
1793	91	if (present(dim_names)) vd%z_grid = '1' ! In this case the names are used to set the non-horizontal axes
1794	91	if (present(fixed)) then ; if (fixed) vd%t_grid = '1' ; endif
1795
1796	91	vd%cmor_field_name = ""
1797	91	vd%cmor_units = ""
1798	91	vd%cmor_longname = ""
1799	91	vd%conversion = 1.0
1800	546	vd%dim_names(:) = ""
1801
1802		call modify_vardesc(vd, units=units, longname=longname, hor_grid=hor_grid, &
1803		z_grid=z_grid, t_grid=t_grid, position=position, dim_names=dim_names, &
1804		cmor_field_name=cmor_field_name, cmor_units=cmor_units, &
1805		cmor_longname=cmor_longname, conversion=conversion, caller=cllr, &
1806	91	extra_axes=extra_axes)
1807
1808	637	end function var_desc
1809
1810
1811		!> This routine modifies the named elements of a vardesc type.
1812		!! All arguments are optional, except the vardesc type to be modified.
1813	93	subroutine modify_vardesc(vd, name, units, longname, hor_grid, z_grid, t_grid, &
1814	93	cmor_field_name, cmor_units, cmor_longname, conversion, caller, position, dim_names, &
1815	93	extra_axes)
1816		type(vardesc), intent(inout) :: vd !< vardesc type that is modified
1817		character(len=*), optional, intent(in) :: name !< name of variable
1818		character(len=*), optional, intent(in) :: units !< units of variable
1819		character(len=*), optional, intent(in) :: longname !< long name of variable
1820		character(len=*), optional, intent(in) :: hor_grid !< horizontal staggering of variable
1821		character(len=*), optional, intent(in) :: z_grid !< vertical staggering of variable
1822		character(len=*), optional, intent(in) :: t_grid !< time description: s, p, or 1
1823		character(len=*), optional, intent(in) :: cmor_field_name !< CMOR name
1824		character(len=*), optional, intent(in) :: cmor_units !< CMOR physical dimensions of variable
1825		character(len=*), optional, intent(in) :: cmor_longname !< CMOR long name
1826		real , optional, intent(in) :: conversion !< A multiplicative factor for unit conversions,
1827		!! such as needed to convert from intensive to
1828		!! extensive or dimensional consistency testing
1829		!! [various] or [a A-1 ~> 1]
1830		character(len=*), optional, intent(in) :: caller !< The calling routine for error messages
1831		integer, optional, intent(in) :: position !< A coded integer indicating the horizontal position
1832		!! of this variable if it has such dimensions.
1833		!! Valid values include CORNER, CENTER, EAST_FACE
1834		!! NORTH_FACE, and 0 for no horizontal dimensions.
1835		character(len=*), dimension(:), &
1836		optional, intent(in) :: dim_names !< The names of the dimensions of this variable
1837		type(axis_info), dimension(:), &
1838		optional, intent(in) :: extra_axes !< dimensions other than space-time
1839
1840		character(len=120) :: cllr
1841		integer :: n
1842
1843	93	cllr = "mod_vardesc" ; if (present(caller)) cllr = trim(caller)
1844
1845	93	if (present(name)) call safe_string_copy(name, vd%name, "vd%name", cllr)
1846
1847	93	if (present(longname)) call safe_string_copy(longname, vd%longname, &
1848	91	"vd%longname of "//trim(vd%name), cllr)
1849	93	if (present(units)) call safe_string_copy(units, vd%units, &
1850	91	"vd%units of "//trim(vd%name), cllr)
1851	93	if (present(position)) then
1852	0	vd%position = position
1853	0	select case (position)
1854	0	case (CENTER) ; vd%hor_grid = 'T'
1855	0	case (CORNER) ; vd%hor_grid = 'Bu'
1856	0	case (EAST_FACE) ; vd%hor_grid = 'Cu'
1857	0	case (NORTH_FACE) ; vd%hor_grid = 'Cv'
1858	0	case (0) ; vd%hor_grid = '1'
1859		case default
1860	0	call MOM_error(FATAL, "modify_vardesc: "//trim(vd%name)//" has unrecognized position argument")
1861		end select
1862		endif
1863	93	if (present(hor_grid)) then
1864	69	call safe_string_copy(hor_grid, vd%hor_grid, "vd%hor_grid of "//trim(vd%name), cllr)
1865	69	vd%position = position_from_horgrid(vd%hor_grid)
1866	69	if (present(caller) .and. (vd%position == -1)) then
1867		call MOM_error(FATAL, "modify_vardesc called by "//trim(caller)//": "//trim(vd%name)//&
1868	0	" has an unrecognized hor_grid argument "//trim(vd%hor_grid))
1869	69	elseif (vd%position == -1) then
1870	0	call MOM_error(FATAL, "modify_vardesc called with bad hor_grid argument "//trim(vd%hor_grid))
1871		endif
1872		endif
1873	93	if (present(z_grid)) call safe_string_copy(z_grid, vd%z_grid, &
1874	82	"vd%z_grid of "//trim(vd%name), cllr)
1875	93	if (present(t_grid)) call safe_string_copy(t_grid, vd%t_grid, &
1876	4	"vd%t_grid of "//trim(vd%name), cllr)
1877
1878	93	if (present(cmor_field_name)) call safe_string_copy(cmor_field_name, vd%cmor_field_name, &
1879	3	"vd%cmor_field_name of "//trim(vd%name), cllr)
1880	93	if (present(cmor_units)) call safe_string_copy(cmor_units, vd%cmor_units, &
1881	0	"vd%cmor_units of "//trim(vd%name), cllr)
1882	93	if (present(cmor_longname)) call safe_string_copy(cmor_longname, vd%cmor_longname, &
1883	2	"vd%cmor_longname of "//trim(vd%name), cllr)
1884
1885	93	if (present(conversion)) vd%conversion = conversion
1886
1887	93	if (present(dim_names)) then
1888	0	do n=1,min(5,size(dim_names)) ; if (len_trim(dim_names(n)) > 0) then
1889	0	call safe_string_copy(dim_names(n), vd%dim_names(n), "vd%dim_names of "//trim(vd%name), cllr)
1890		endif ; enddo
1891		endif
1892
1893	93	if (present(extra_axes)) then
1894	0	do n=1,size(extra_axes) ; if (len_trim(extra_axes(n)%name) > 0) then
1895	0	vd%extra_axes(n) = extra_axes(n)
1896		endif ; enddo
1897		endif
1898
1899	93	end subroutine modify_vardesc
1900
1901	69	integer function position_from_horgrid(hor_grid)
1902		character(len=*), intent(in) :: hor_grid !< horizontal staggering of variable
1903
1904	69	select case (trim(hor_grid))
1905	0	case ('h') ; position_from_horgrid = CENTER
1906	0	case ('q') ; position_from_horgrid = CORNER
1907	4	case ('u') ; position_from_horgrid = EAST_FACE
1908	4	case ('v') ; position_from_horgrid = NORTH_FACE
1909	0	case ('T') ; position_from_horgrid = CENTER
1910	1	case ('Bu') ; position_from_horgrid = CORNER
1911	1	case ('Cu') ; position_from_horgrid = EAST_FACE
1912	1	case ('Cv') ; position_from_horgrid = NORTH_FACE
1913	58	case ('1') ; position_from_horgrid = 0
1914	138	case default ; position_from_horgrid = -1 ! This is a bad-value flag.
1915		end select
1916	138	end function position_from_horgrid
1917
1918		!> Store information that can be used to create an axis in a subsequent call to create_file.
1919	315	subroutine set_axis_info(axis, name, units, longname, ax_size, ax_data, cartesian, sense)
1920		type(axis_info), intent(inout) :: axis !< A type with information about a named axis
1921		character(len=*), intent(in) :: name !< The name of this axis for use in files
1922		character(len=*), optional, intent(in) :: units !< The units of the axis labels
1923		character(len=*), optional, intent(in) :: longname !< Long name of the axis variable
1924		integer, optional, intent(in) :: ax_size !< The number of elements in this axis
1925		real, dimension(:), optional, intent(in) :: ax_data !< The values of the data on the axis [arbitrary]
1926		character(len=*), optional, intent(in) :: cartesian !< A variable indicating which direction this axis
1927		!! axis corresponds with. Valid values
1928		!! include 'X', 'Y', 'Z', 'T', and 'N' (the default) for none.
1929		integer, optional, intent(in) :: sense !< This is 1 for axes whose values increase upward, or -1
1930		!! if they increase downward. The default, 0, is ignored.
1931
1932	315	call safe_string_copy(name, axis%name, "axis%name of "//trim(name), "set_axis_info")
1933		! Set the default values.
1934	315	axis%longname = trim(axis%name) ; axis%units = "" ; axis%cartesian = "N" ; axis%sense = 0
1935
1936	315	if (present(longname)) call safe_string_copy(longname, axis%longname, &
1937	0	"axis%longname of "//trim(name), "set_axis_info")
1938	315	if (present(units)) call safe_string_copy(units, axis%units, &
1939	0	"axis%units of "//trim(name), "set_axis_info")
1940	315	if (present(cartesian)) call safe_string_copy(cartesian, axis%cartesian, &
1941	0	"axis%cartesian of "//trim(name), "set_axis_info")
1942	315	if (present(sense)) axis%sense = sense
1943
1944	315	if (.not.(present(ax_size) .or. present(ax_data)) ) then
1945		call MOM_error(FATAL, "set_axis_info called for "//trim(name)//&
1946	0	"without either an ax_size or an ax_data argument.")
1947	315	elseif (present(ax_size) .and. present(ax_data)) then
1948	0	if (size(ax_data) /= ax_size) call MOM_error(FATAL, "set_axis_info called for "//trim(name)//&
1949	0	"with an inconsistent value of ax_size and size of ax_data.")
1950		endif
1951
1952	315	if (present(ax_size)) then
1953	315	axis%ax_size = ax_size
1954		else
1955	0	axis%ax_size = size(ax_data)
1956		endif
1957	315	if (present(ax_data)) then
1958	0	allocate(axis%ax_data(axis%ax_size)) ; axis%ax_data(:) = ax_data(:)
1959		endif
1960
1961	315	end subroutine set_axis_info
1962
1963		!> Delete the information in an array of axis_info types and deallocate memory in them.
1964	0	subroutine delete_axis_info(axes)
1965		type(axis_info), dimension(:), intent(inout) :: axes !< An array with information about named axes
1966
1967		integer :: n
1968	0	do n=1,size(axes)
1969	0	axes(n)%name = "" ; axes(n)%longname = "" ; axes(n)%units = "" ; axes(n)%cartesian = "N"
1970	0	axes(n)%sense = 0 ; axes(n)%ax_size = 0
1971	0	if (allocated(axes(n)%ax_data)) deallocate(axes(n)%ax_data)
1972		enddo
1973	0	end subroutine delete_axis_info
1974
1975
1976		!> Retrieve the information from an axis_info type.
1977	0	subroutine get_axis_info(axis,name,longname,units,cartesian,ax_size,ax_data)
1978		type(axis_info), intent(in) :: axis !< An axis type
1979		character(len=*), intent(out), optional :: name !< The axis name.
1980		character(len=*), intent(out), optional :: longname !< The axis longname.
1981		character(len=*), intent(out), optional :: units !< The axis units.
1982		character(len=*), intent(out), optional :: cartesian !< The cartesian attribute
1983		!! of the axis [X,Y,Z,T].
1984		integer, intent(out), optional :: ax_size !< The size of the axis.
1985		real, optional, allocatable, dimension(:), intent(out) :: ax_data !< The axis label data [arbitrary]
1986
1987	0	if (present(ax_data)) then
1988	0	if (allocated(ax_data)) deallocate(ax_data)
1989	0	allocate(ax_data(axis%ax_size))
1990	0	ax_data(:) = axis%ax_data
1991		endif
1992
1993	0	if (present(name)) name = axis%name
1994	0	if (present(longname)) longname = axis%longname
1995	0	if (present(units)) units = axis%units
1996	0	if (present(cartesian)) cartesian = axis%cartesian
1997	0	if (present(ax_size)) ax_size = axis%ax_size
1998
1999	0	end subroutine get_axis_info
2000
2001		!> Store information that can be used to create an attribute in a subsequent call to create_file.
2002	0	subroutine set_attribute_info(attribute, name, str_value)
2003		type(attribute_info), intent(inout) :: attribute !< A type with information about a named attribute
2004		character(len=*), intent(in) :: name !< The name of this attribute for use in files
2005		character(len=*), intent(in) :: str_value !< The value of this attribute
2006
2007	0	attribute%name = trim(name)
2008	0	attribute%att_val = trim(str_value)
2009	0	end subroutine set_attribute_info
2010
2011		!> Delete the information in an array of attribute_info types and deallocate memory in them.
2012	0	subroutine delete_attribute_info(atts)
2013		type(attribute_info), dimension(:), intent(inout) :: atts !< An array of global attributes
2014
2015		integer :: n
2016	0	do n=1,size(atts)
2017	0	if (allocated(atts(n)%name)) deallocate(atts(n)%name)
2018	0	if (allocated(atts(n)%att_val)) deallocate(atts(n)%att_val)
2019		enddo
2020	0	end subroutine delete_attribute_info
2021
2022
2023		!> This function returns the CMOR standard name given a CMOR longname, based on
2024		!! the standard pattern of character conversions.
2025	11	function cmor_long_std(longname) result(std_name)
2026		character(len=*), intent(in) :: longname !< The CMOR longname being converted
2027		character(len=len(longname)) :: std_name !< The CMOR standard name generated from longname
2028
2029		integer :: k
2030
2031	11	std_name = lowercase(longname)
2032
2033	798	do k=1, len_trim(std_name)
2034	798	if (std_name(k:k) == ' ') std_name(k:k) = '_'
2035		enddo
2036
2037	11	end function cmor_long_std
2038
2039		!> This routine queries vardesc
2040	861	subroutine query_vardesc(vd, name, units, longname, hor_grid, z_grid, t_grid, &
2041		cmor_field_name, cmor_units, cmor_longname, conversion, caller, &
2042	168	extra_axes, position, dim_names)
2043		type(vardesc), intent(in) :: vd !< vardesc type that is queried
2044		character(len=*), optional, intent(out) :: name !< name of variable
2045		character(len=*), optional, intent(out) :: units !< units of variable
2046		character(len=*), optional, intent(out) :: longname !< long name of variable
2047		character(len=*), optional, intent(out) :: hor_grid !< horizontal staggering of variable
2048		character(len=*), optional, intent(out) :: z_grid !< verticle staggering of variable
2049		character(len=*), optional, intent(out) :: t_grid !< time description: s, p, or 1
2050		character(len=*), optional, intent(out) :: cmor_field_name !< CMOR name
2051		character(len=*), optional, intent(out) :: cmor_units !< CMOR physical dimensions of variable
2052		character(len=*), optional, intent(out) :: cmor_longname !< CMOR long name
2053		real , optional, intent(out) :: conversion !< for unit conversions, such as needed to
2054		!! convert from intensive to extensive
2055		!! [various] or [a A-1 ~> 1]
2056		character(len=*), optional, intent(in) :: caller !< calling routine?
2057		type(axis_info), dimension(5), &
2058		optional, intent(out) :: extra_axes !< dimensions other than space-time
2059		integer, optional, intent(out) :: position !< A coded integer indicating the horizontal position
2060		!! of this variable if it has such dimensions.
2061		!! Valid values include CORNER, CENTER, EAST_FACE
2062		!! NORTH_FACE, and 0 for no horizontal dimensions.
2063		character(len=*), dimension(:), &
2064		optional, intent(out) :: dim_names !< The names of the dimensions of this variable
2065
2066		integer :: n
2067		integer, parameter :: nmax_extraaxes = 5
2068		character(len=120) :: cllr, varname
2069	168	cllr = "mod_vardesc"
2070	168	if (present(caller)) cllr = trim(caller)
2071
2072	168	if (present(name)) call safe_string_copy(vd%name, name, &
2073	103	"vd%name of "//trim(vd%name), cllr)
2074	168	if (present(longname)) call safe_string_copy(vd%longname, longname, &
2075	3	"vd%longname of "//trim(vd%name), cllr)
2076	168	if (present(units)) call safe_string_copy(vd%units, units, &
2077	6	"vd%units of "//trim(vd%name), cllr)
2078	168	if (present(hor_grid)) call safe_string_copy(vd%hor_grid, hor_grid, &
2079	0	"vd%hor_grid of "//trim(vd%name), cllr)
2080	168	if (present(z_grid)) call safe_string_copy(vd%z_grid, z_grid, &
2081	63	"vd%z_grid of "//trim(vd%name), cllr)
2082	168	if (present(t_grid)) call safe_string_copy(vd%t_grid, t_grid, &
2083	65	"vd%t_grid of "//trim(vd%name), cllr)
2084
2085	168	if (present(cmor_field_name)) call safe_string_copy(vd%cmor_field_name, cmor_field_name, &
2086	3	"vd%cmor_field_name of "//trim(vd%name), cllr)
2087	168	if (present(cmor_units)) call safe_string_copy(vd%cmor_units, cmor_units, &
2088	0	"vd%cmor_units of "//trim(vd%name), cllr)
2089	168	if (present(cmor_longname)) call safe_string_copy(vd%cmor_longname, cmor_longname, &
2090	3	"vd%cmor_longname of "//trim(vd%name), cllr)
2091
2092	168	if (present(conversion)) conversion = vd%conversion
2093
2094	168	if (present(position)) then
2095	128	position = vd%position
2096	128	if (position == -1) position = position_from_horgrid(vd%hor_grid)
2097		endif
2098	168	if (present(dim_names)) then
2099	0	do n=1,min(5,size(dim_names))
2100	0	call safe_string_copy(vd%dim_names(n), dim_names(n), "vd%dim_names of "//trim(vd%name), cllr)
2101		enddo
2102		endif
2103
2104	168	if (present(extra_axes)) then
2105		! save_restart expects 5 extra axes (can be empty)
2106	378	do n=1, nmax_extraaxes
2107	378	if (vd%extra_axes(n)%ax_size>=1) then
2108	0	extra_axes(n) = vd%extra_axes(n)
2109		else
2110		! return an empty axis
2111	315	write(varname,"('dummy',i1.1)") n
2112	315	call set_axis_info(extra_axes(n), name=trim(varname), ax_size=1)
2113		endif
2114		enddo
2115		endif
2116
2117	168	end subroutine query_vardesc
2118
2119
2120		!> Read a scalar from file using infrastructure I/O.
2121	0	subroutine MOM_read_data_0d(filename, fieldname, data, timelevel, scale, MOM_Domain, &
2122		global_file, file_may_be_4d)
2123		character(len=*), intent(in) :: filename !< Input filename
2124		character(len=*), intent(in) :: fieldname !< Field variable name
2125		real, intent(inout) :: data !< Field value in arbitrary units [A ~> a]
2126		integer, optional, intent(in) :: timelevel !< Time level to read in file
2127		real, optional, intent(in) :: scale !< A scaling factor that the variable is multiplied by
2128		!! before it is returned to convert from the units in the file
2129		!! to the internal units for this variable [A a-1 ~> 1]
2130		type(MOM_domain_type), optional, intent(in) :: MOM_Domain !< Model domain decomposition
2131		logical, optional, intent(in) :: global_file !< If true, read from a single file
2132		logical, optional, intent(in) :: file_may_be_4d !< If true, fields may be stored
2133		!! as 4d arrays in the file.
2134
2135		call read_field(filename, fieldname, data, &
2136		timelevel=timelevel, scale=scale, MOM_Domain=MOM_Domain, &
2137	0	global_file=global_file, file_may_be_4d=file_may_be_4d)
2138	0	end subroutine MOM_read_data_0d
2139
2140
2141		!> Read a scalar integer from file using infrastructure I/O.
2142	0	subroutine MOM_read_data_0d_int(filename, fieldname, data, timelevel)
2143		character(len=*), intent(in) :: filename !< Input filename
2144		character(len=*), intent(in) :: fieldname !< Field variable name
2145		integer, intent(inout) :: data !< Field value
2146		integer, optional, intent(in) :: timelevel !< Time level to read in file
2147
2148	0	call read_field(filename, fieldname, data, timelevel=timelevel)
2149	0	end subroutine MOM_read_data_0d_int
2150
2151
2152		!> Read a 1d array from file using infrastructure I/O.
2153	0	subroutine MOM_read_data_1d(filename, fieldname, data, timelevel, scale, MOM_Domain, &
2154		global_file, file_may_be_4d)
2155		character(len=*), intent(in) :: filename !< Input filename
2156		character(len=*), intent(in) :: fieldname !< Field variable name
2157		real, dimension(:), intent(inout) :: data !< Field value in arbitrary units [A ~> a]
2158		integer, optional, intent(in) :: timelevel !< Time level to read in file
2159		real, optional, intent(in) :: scale !< A scaling factor that the variable is multiplied by
2160		!! before it is returned to convert from the units in the file
2161		!! to the internal units for this variable [A a-1 ~> 1]
2162		type(MOM_domain_type), optional, intent(in) :: MOM_Domain !< Model domain decomposition
2163		logical, optional, intent(in) :: global_file !< If true, read from a single file
2164		logical, optional, intent(in) :: file_may_be_4d !< If true, fields may be stored
2165		!! as 4d arrays in the file.
2166
2167		call read_field(filename, fieldname, data, &
2168		timelevel=timelevel, scale=scale, MOM_Domain=MOM_Domain, &
2169	0	global_file=global_file, file_may_be_4d=file_may_be_4d)
2170
2171	0	end subroutine MOM_read_data_1d
2172
2173
2174		!> Read a 1d integer array from file using infrastructure I/O.
2175	0	subroutine MOM_read_data_1d_int(filename, fieldname, data, timelevel)
2176		character(len=*), intent(in) :: filename !< Input filename
2177		character(len=*), intent(in) :: fieldname !< Field variable name
2178		integer, dimension(:), intent(inout) :: data !< Field value
2179		integer, optional, intent(in) :: timelevel !< Time level to read in file
2180
2181	0	call read_field(filename, fieldname, data, timelevel=timelevel)
2182	0	end subroutine MOM_read_data_1d_int
2183
2184
2185		!> Read a 2d array from file using infrastructure I/O.
2186	0	subroutine MOM_read_data_2d(filename, fieldname, data, MOM_Domain, timelevel, position, &
2187		scale, global_file, file_may_be_4d, turns)
2188		character(len=*), intent(in) :: filename !< Input filename
2189		character(len=*), intent(in) :: fieldname !< Field variable name
2190		real, dimension(:,:), intent(inout) :: data !< Field value in arbitrary units [A ~> a]
2191		type(MOM_domain_type), target, &
2192		intent(in) :: MOM_Domain !< Model domain decomposition
2193		integer, optional, intent(in) :: timelevel !< Time level to read in file
2194		integer, optional, intent(in) :: position !< Grid positioning flag
2195		real, optional, intent(in) :: scale !< A scaling factor that the variable is multiplied by
2196		!! before it is returned to convert from the units in the file
2197		!! to the internal units for this variable [A a-1 ~> 1]
2198		logical, optional, intent(in) :: global_file !< If true, read from a single file
2199		logical, optional, intent(in) :: file_may_be_4d !< If true, fields may be stored
2200		!! as 4d arrays in the file.
2201		integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data. If absent
2202		!! the number of turns is taken from MOM_Domain.
2203
2204		! Local variables
2205		integer :: qturns ! Number of quarter-turns from input to model grid
2206	0	real, allocatable :: data_in(:,:) ! Field array on the input grid in arbitrary units [A ~> a]
2207		type(MOM_domain_type), pointer :: domain_ptr => NULL() ! Pointer to the unrotated domain for reading
2208
2209	0	qturns = MOM_domain%turns ; if (present(turns)) qturns = modulo(turns, 4)
2210
2211	0	domain_ptr => MOM_Domain
2212	0	if (associated(MOM_Domain%domain_in) .and. (qturns /= 0)) domain_ptr => MOM_Domain%domain_in
2213
2214	0	if (qturns == 0) then
2215		call read_field(filename, fieldname, data, MOM_Domain, &
2216		timelevel=timelevel, position=position, scale=scale, &
2217	0	global_file=global_file, file_may_be_4d=file_may_be_4d)
2218		else
2219	0	call allocate_rotated_array(data, [1,1], -qturns, data_in)
2220	0	call rotate_array(data, -qturns, data_in)
2221		call read_field(filename, fieldname, data_in, domain_ptr, &
2222		timelevel=timelevel, position=position, scale=scale, &
2223	0	global_file=global_file, file_may_be_4d=file_may_be_4d)
2224	0	call rotate_array(data_in, qturns, data)
2225	0	deallocate(data_in)
2226		endif
2227
2228	0	end subroutine MOM_read_data_2d
2229
2230
2231		!> Read a 2d array (which might have halos) from a file using native netCDF I/O.
2232	0	subroutine read_netCDF_data_2d(filename, fieldname, values, MOM_Domain, &
2233		timelevel, position, rescale, turns)
2234		character(len=*), intent(in) :: filename
2235		!< Input filename
2236		character(len=*), intent(in) :: fieldname
2237		!< Field variable name
2238		real, intent(inout) :: values(:,:)
2239		!< Field values read from the file. It would be intent(out) but for the
2240		!! need to preserve any initialized values in the halo regions.
2241		type(MOM_domain_type), intent(in) :: MOM_Domain
2242		!< Model domain decomposition
2243		integer, optional, intent(in) :: timelevel
2244		!< Time level to read in file
2245		integer, optional, intent(in) :: position
2246		!< Grid positioning flag
2247		real, optional, intent(in) :: rescale
2248		!< Rescale factor, omitting this is the same as setting it to 1.
2249		integer, optional, intent(in) :: turns
2250		!< Number of quarter-turns to rotate the data. If absent the number of turns is taken
2251		!! from MOM_Domain.
2252
2253		integer :: qturns
2254		! Number of quarter-turns from input to model grid
2255	0	real, allocatable :: values_in(:,:)
2256		! Field array on the unrotated input grid
2257	0	type(MOM_netcdf_file) :: handle
2258		! netCDF file handle
2259
2260		! General-purpose IO will require the following arguments, but they are not
2261		! yet implemented, so we raise an error if they are present.
2262
2263		! Fields are currently assumed on cell centers, and position is unsupported
2264	0	if (present(position)) &
2265	0	call MOM_error(FATAL, 'read_netCDF_data: position is not yet supported.')
2266
2267		! Timelevels are not yet supported
2268	0	if (present(timelevel)) &
2269	0	call MOM_error(FATAL, 'read_netCDF_data: timelevel is not yet supported.')
2270
2271	0	call handle%open(filename, action=READONLY_FILE, MOM_domain=MOM_domain)
2272	0	call handle%update()
2273
2274	0	qturns = MOM_domain%turns ; if (present(turns)) qturns = modulo(turns, 4)
2275
2276	0	if (qturns == 0) then
2277	0	call handle%read(fieldname, values, rescale=rescale)
2278		else
2279	0	call allocate_rotated_array(values, [1,1], -qturns, values_in)
2280	0	call rotate_array(values, -qturns, values_in)
2281	0	call handle%read(fieldname, values_in, rescale=rescale)
2282	0	call rotate_array(values_in, qturns, values)
2283	0	deallocate(values_in)
2284		endif
2285
2286	0	call handle%close()
2287	0	end subroutine read_netCDF_data_2d
2288
2289
2290		!> Read a 2d region array from file using infrastructure I/O.
2291	0	subroutine MOM_read_data_2d_region(filename, fieldname, data, start, nread, MOM_domain, &
2292		no_domain, scale, turns)
2293		character(len=*), intent(in) :: filename !< Input filename
2294		character(len=*), intent(in) :: fieldname !< Field variable name
2295		real, dimension(:,:), intent(inout) :: data !< Field value in arbitrary units [A ~> a]
2296		integer, dimension(:), intent(in) :: start !< Starting index for each axis.
2297		!! In 2d, start(3:4) must be 1.
2298		integer, dimension(:), intent(in) :: nread !< Number of values to read along each axis.
2299		!! In 2d, nread(3:4) must be 1.
2300		type(MOM_domain_type), optional, intent(in) :: MOM_Domain !< Model domain decomposition
2301		logical, optional, intent(in) :: no_domain !< If true, field does not use
2302		!! domain decomposion.
2303		real, optional, intent(in) :: scale !< A scaling factor that the variable is multiplied by
2304		!! before it is returned to convert from the units in the file
2305		!! to the internal units for this variable [A a-1 ~> 1]
2306		integer, optional, intent(in) :: turns !< Number of quarter turns from
2307		!! input to model grid
2308
2309		integer :: qturns ! Number of quarter turns
2310	0	real, allocatable :: data_in(:,:) ! Field array on the input grid in arbitrary units [A ~> a]
2311
2312	0	qturns = 0
2313	0	if (present(turns)) qturns = modulo(turns, 4)
2314
2315	0	if (qturns == 0) then
2316		call read_field(filename, fieldname, data, start, nread, &
2317	0	MOM_Domain=MOM_Domain, no_domain=no_domain, scale=scale)
2318		else
2319	0	call allocate_rotated_array(data, [1,1], -qturns, data_in)
2320	0	call rotate_array(data, -qturns, data_in)
2321	0	if (associated(MOM_Domain%domain_in)) then
2322		call read_field(filename, fieldname, data_in, start, nread, &
2323	0	MOM_Domain=MOM_Domain%domain_in, no_domain=no_domain, scale=scale)
2324		else
2325		call read_field(filename, fieldname, data_in, start, nread, &
2326	0	MOM_Domain=MOM_Domain, no_domain=no_domain, scale=scale)
2327		endif
2328	0	call rotate_array(data_in, qturns, data)
2329	0	deallocate(data_in)
2330		endif
2331	0	end subroutine MOM_read_data_2d_region
2332
2333
2334		!> Read a 3d array from file using infrastructure I/O.
2335	0	subroutine MOM_read_data_3d(filename, fieldname, data, MOM_Domain, timelevel, position, &
2336		scale, global_file, file_may_be_4d, turns)
2337		character(len=*), intent(in) :: filename !< Input filename
2338		character(len=*), intent(in) :: fieldname !< Field variable name
2339		real, dimension(:,:,:), intent(inout) :: data !< Field value in arbitrary units [A ~> a]
2340		type(MOM_domain_type), target, &
2341		intent(in) :: MOM_Domain !< Model domain decomposition
2342		integer, optional, intent(in) :: timelevel !< Time level to read in file
2343		integer, optional, intent(in) :: position !< Grid positioning flag
2344		real, optional, intent(in) :: scale !< A scaling factor that the variable is multiplied by
2345		!! before it is returned to convert from the units in the file
2346		!! to the internal units for this variable [A a-1 ~> 1]
2347		logical, optional, intent(in) :: global_file !< If true, read from a single file
2348		logical, optional, intent(in) :: file_may_be_4d !< If true, fields may be stored
2349		!! as 4d arrays in the file.
2350		integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data. If absent
2351		!! the number of turns is taken from MOM_Domain.
2352
2353		! Local variables
2354		integer :: qturns ! Number of quarter-turns from input to model grid
2355	0	real, allocatable :: data_in(:,:,:) ! Field array on the input grid in arbitrary units [A ~> a]
2356		type(MOM_domain_type), pointer :: domain_ptr => NULL() ! Pointer to the unrotated domain for reading
2357
2358	0	domain_ptr => MOM_Domain
2359	0	if (associated(MOM_Domain%domain_in) .and. (qturns /= 0)) domain_ptr => MOM_Domain%domain_in
2360
2361	0	qturns = MOM_domain%turns ; if (present(turns)) qturns = modulo(turns, 4)
2362	0	if (qturns == 0) then
2363		call read_field(filename, fieldname, data, MOM_Domain, &
2364		timelevel=timelevel, position=position, scale=scale, &
2365	0	global_file=global_file, file_may_be_4d=file_may_be_4d)
2366		else
2367	0	call allocate_rotated_array(data, [1,1,1], -qturns, data_in)
2368	0	call rotate_array(data, -qturns, data_in)
2369		call read_field(filename, fieldname, data_in, domain_ptr, &
2370		timelevel=timelevel, position=position, scale=scale, &
2371	0	global_file=global_file, file_may_be_4d=file_may_be_4d)
2372	0	call rotate_array(data_in, qturns, data)
2373	0	deallocate(data_in)
2374		endif
2375
2376	0	end subroutine MOM_read_data_3d
2377
2378		!> Read a 3d region array from file using infrastructure I/O.
2379	0	subroutine MOM_read_data_3d_region(filename, fieldname, data, start, nread, MOM_domain, &
2380		no_domain, scale, turns)
2381		character(len=*), intent(in) :: filename !< Input filename
2382		character(len=*), intent(in) :: fieldname !< Field variable name
2383		real, dimension(:,:,:), intent(inout) :: data !< Field value in arbitrary units [A ~> a]
2384		integer, dimension(:), intent(in) :: start !< Starting index for each axis.
2385		integer, dimension(:), intent(in) :: nread !< Number of values to read along each axis.
2386		type(MOM_domain_type), optional, intent(in) :: MOM_Domain !< Model domain decomposition
2387		logical, optional, intent(in) :: no_domain !< If true, field does not use
2388		!! domain decomposion.
2389		real, optional, intent(in) :: scale !< A scaling factor that the variable is multiplied by
2390		!! before it is returned to convert from the units in the file
2391		!! to the internal units for this variable [A a-1 ~> 1]
2392		integer, optional, intent(in) :: turns !< Number of quarter turns from
2393		!! input to model grid
2394
2395		integer :: qturns ! Number of quarter turns
2396	0	real, allocatable :: data_in(:,:,:) ! Field array on the input grid in arbitrary units [A ~> a]
2397
2398	0	qturns = 0
2399	0	if (present(turns)) qturns = modulo(turns, 4)
2400
2401	0	if (qturns == 0) then
2402		call read_field(filename, fieldname, data, start, nread, &
2403	0	MOM_Domain=MOM_Domain, no_domain=no_domain, scale=scale)
2404		else
2405	0	call allocate_rotated_array(data, [1,1,1], -qturns, data_in)
2406	0	call rotate_array(data, -qturns, data_in)
2407	0	if (associated(MOM_Domain%domain_in)) then
2408		call read_field(filename, fieldname, data_in, start, nread, &
2409	0	MOM_Domain=MOM_Domain%domain_in, no_domain=no_domain, scale=scale)
2410		else
2411		call read_field(filename, fieldname, data_in, start, nread, &
2412	0	MOM_Domain=MOM_Domain, no_domain=no_domain, scale=scale)
2413		endif
2414	0	call rotate_array(data_in, qturns, data)
2415	0	deallocate(data_in)
2416		endif
2417	0	end subroutine MOM_read_data_3d_region
2418
2419		!> Read a 4d array from file using infrastructure I/O.
2420	0	subroutine MOM_read_data_4d(filename, fieldname, data, MOM_Domain, &
2421		timelevel, position, scale, global_file, turns)
2422		character(len=*), intent(in) :: filename !< Input filename
2423		character(len=*), intent(in) :: fieldname !< Field variable name
2424		real, dimension(:,:,:,:), intent(inout) :: data !< Field value in arbitrary units [A ~> a]
2425		type(MOM_domain_type), target, &
2426		intent(in) :: MOM_Domain !< Model domain decomposition
2427		integer, optional, intent(in) :: timelevel !< Time level to read in file
2428		integer, optional, intent(in) :: position !< Grid positioning flag
2429		real, optional, intent(in) :: scale !< A scaling factor that the variable is multiplied by
2430		!! before it is returned to convert from the units in the file
2431		!! to the internal units for this variable [A a-1 ~> 1]
2432		logical, optional, intent(in) :: global_file !< If true, read from a single file
2433		integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data. If absent
2434		!! the number of turns is taken from MOM_Domain.
2435
2436		! Local variables
2437		integer :: qturns ! Number of quarter-turns from input to model grid
2438	0	real, allocatable :: data_in(:,:,:,:) ! Field array on the input grid in arbitrary units [A ~> a]
2439		type(MOM_domain_type), pointer :: domain_ptr => NULL() ! Pointer to the unrotated domain for reading
2440
2441	0	qturns = MOM_domain%turns ; if (present(turns)) qturns = modulo(turns, 4)
2442
2443	0	domain_ptr => MOM_Domain
2444	0	if (associated(MOM_Domain%domain_in) .and. (qturns /= 0)) domain_ptr => MOM_Domain%domain_in
2445
2446	0	if (qturns == 0) then
2447		call read_field(filename, fieldname, data, MOM_Domain, &
2448		timelevel=timelevel, position=position, scale=scale, &
2449	0	global_file=global_file)
2450		else
2451		! Read field along the input grid and rotate to the model grid
2452	0	call allocate_rotated_array(data, [1,1,1,1], -qturns, data_in)
2453	0	call rotate_array(data, -qturns, data_in)
2454		call read_field(filename, fieldname, data_in, domain_ptr, timelevel=timelevel, &
2455	0	position=position, scale=scale, global_file=global_file)
2456	0	call rotate_array(data_in, qturns, data)
2457	0	deallocate(data_in)
2458		endif
2459
2460	0	end subroutine MOM_read_data_4d
2461
2462
2463		!> Read a 2d vector tuple from file using infrastructure I/O.
2464	0	subroutine MOM_read_vector_2d(filename, u_fieldname, v_fieldname, u_data, v_data, MOM_Domain, &
2465		timelevel, stagger, scalar_pair, scale, turns)
2466		character(len=*), intent(in) :: filename !< Input filename
2467		character(len=*), intent(in) :: u_fieldname !< Field variable name in u
2468		character(len=*), intent(in) :: v_fieldname !< Field variable name in v
2469		real, dimension(:,:), intent(inout) :: u_data !< Field value at u points in arbitrary units [A ~> a]
2470		real, dimension(:,:), intent(inout) :: v_data !< Field value at v points in arbitrary units [A ~> a]
2471		type(MOM_domain_type), target, &
2472		intent(in) :: MOM_Domain !< Model domain decomposition
2473		integer, optional, intent(in) :: timelevel !< Time level to read in file
2474		integer, optional, intent(in) :: stagger !< Grid staggering flag
2475		logical, optional, intent(in) :: scalar_pair !< True if tuple is not a vector
2476		real, optional, intent(in) :: scale !< A scaling factor that the vector is multiplied by
2477		!! before it is returned to convert from the units in the file
2478		!! to the internal units for this variable [A a-1 ~> 1]
2479		integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data. If absent
2480		!! the number of turns is taken from MOM_Domain.
2481
2482		! Local variables
2483		integer :: qturns ! Number of quarter-turns from input to model grid
2484	0	real, allocatable :: u_data_in(:,:), v_data_in(:,:) ! [uv] on the input grid in arbitrary units [A ~> a]
2485		type(MOM_domain_type), pointer :: domain_ptr => NULL() ! Pointer to the unrotated domain for reading
2486
2487	0	qturns = MOM_domain%turns ; if (present(turns)) qturns = modulo(turns, 4)
2488
2489	0	domain_ptr => MOM_Domain
2490	0	if (associated(MOM_Domain%domain_in) .and. (qturns /= 0)) domain_ptr => MOM_Domain%domain_in
2491
2492	0	if (qturns == 0) then
2493		call read_vector(filename, u_fieldname, v_fieldname, &
2494		u_data, v_data, MOM_domain, timelevel=timelevel, stagger=stagger, &
2495	0	scalar_pair=scalar_pair, scale=scale)
2496		else
2497	0	call allocate_rotated_array(u_data, [1,1], -qturns, u_data_in)
2498	0	call allocate_rotated_array(v_data, [1,1], -qturns, v_data_in)
2499	0	if (scalar_pair) then
2500	0	call rotate_array_pair(u_data, v_data, -qturns, u_data_in, v_data_in)
2501		else
2502	0	call rotate_vector(u_data, v_data, -qturns, u_data_in, v_data_in)
2503		endif
2504		call read_vector(filename, u_fieldname, v_fieldname, u_data_in, v_data_in, &
2505		domain_ptr, timelevel=timelevel, &
2506	0	stagger=stagger, scalar_pair=scalar_pair, scale=scale)
2507	0	if (scalar_pair) then
2508	0	call rotate_array_pair(u_data_in, v_data_in, qturns, u_data, v_data)
2509		else
2510	0	call rotate_vector(u_data_in, v_data_in, qturns, u_data, v_data)
2511		endif
2512	0	deallocate(v_data_in)
2513	0	deallocate(u_data_in)
2514		endif
2515
2516	0	end subroutine MOM_read_vector_2d
2517
2518
2519		!> Read a 3d vector tuple from file using infrastructure I/O.
2520	0	subroutine MOM_read_vector_3d(filename, u_fieldname, v_fieldname, u_data, v_data, MOM_Domain, &
2521		timelevel, stagger, scalar_pair, scale, turns)
2522		character(len=*), intent(in) :: filename !< Input filename
2523		character(len=*), intent(in) :: u_fieldname !< Field variable name in u
2524		character(len=*), intent(in) :: v_fieldname !< Field variable name in v
2525		real, dimension(:,:,:), intent(inout) :: u_data !< Field value in u in arbitrary units [A ~> a]
2526		real, dimension(:,:,:), intent(inout) :: v_data !< Field value in v in arbitrary units [A ~> a]
2527		type(MOM_domain_type), target, &
2528		intent(in) :: MOM_Domain !< Model domain decomposition
2529		integer, optional, intent(in) :: timelevel !< Time level to read in file
2530		integer, optional, intent(in) :: stagger !< Grid staggering flag
2531		logical, optional, intent(in) :: scalar_pair !< True if tuple is not a vector
2532		real, optional, intent(in) :: scale !< A scaling factor that the vector is multiplied by
2533		!! before it is returned to convert from the units in the file
2534		!! to the internal units for this variable [A a-1 ~> 1]
2535		integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data. If absent
2536		!! the number of turns is taken from MOM_Domain.
2537
2538		! Local variables
2539		integer :: qturns ! Number of quarter-turns from input to model grid
2540	0	real, allocatable :: u_data_in(:,:,:), v_data_in(:,:,:) ! [uv] on the input grid in arbitrary units [A ~> a]
2541		type(MOM_domain_type), pointer :: domain_ptr => NULL() ! Pointer to the unrotated domain for reading
2542
2543	0	qturns = MOM_domain%turns ; if (present(turns)) qturns = modulo(turns, 4)
2544
2545	0	domain_ptr => MOM_Domain
2546	0	if (associated(MOM_Domain%domain_in) .and. (qturns /= 0)) domain_ptr => MOM_Domain%domain_in
2547
2548	0	if (qturns == 0) then
2549		call read_vector(filename, u_fieldname, v_fieldname, &
2550		u_data, v_data, MOM_domain, timelevel=timelevel, stagger=stagger, &
2551	0	scalar_pair=scalar_pair, scale=scale)
2552		else
2553	0	call allocate_rotated_array(u_data, [1,1,1], -qturns, u_data_in)
2554	0	call allocate_rotated_array(v_data, [1,1,1], -qturns, v_data_in)
2555	0	if (scalar_pair) then
2556	0	call rotate_array_pair(u_data, v_data, -qturns, u_data_in, v_data_in)
2557		else
2558	0	call rotate_vector(u_data, v_data, -qturns, u_data_in, v_data_in)
2559		endif
2560		call read_vector(filename, u_fieldname, v_fieldname, u_data_in, v_data_in, &
2561		domain_ptr, timelevel=timelevel, &
2562	0	stagger=stagger, scalar_pair=scalar_pair, scale=scale)
2563	0	if (scalar_pair) then
2564	0	call rotate_array_pair(u_data_in, v_data_in, qturns, u_data, v_data)
2565		else
2566	0	call rotate_vector(u_data_in, v_data_in, qturns, u_data, v_data)
2567		endif
2568	0	deallocate(v_data_in)
2569	0	deallocate(u_data_in)
2570		endif
2571
2572	0	end subroutine MOM_read_vector_3d
2573
2574		!> Write a 4d field to an output file, potentially with rotation
2575	0	subroutine MOM_write_field_legacy_4d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, &
2576		fill_value, turns, scale, unscale, zero_zeros)
2577		type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing
2578		type(fieldtype), intent(in) :: field_md !< Field type with metadata
2579		type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition
2580		real, dimension(:,:,:,:), intent(inout) :: field !< Unrotated field to write in arbitrary units [A ~> a]
2581		real, optional, intent(in) :: tstamp !< Model timestamp, often in [days]
2582		integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1)
2583		real, optional, intent(in) :: fill_value !< Missing data fill value in the units used in the file [a]
2584		integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data
2585		real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied by before
2586		!! it is written [a A-1 ~> 1], for example to convert it
2587		!! from its internal units to the desired units for output
2588		real, optional, intent(in) :: unscale !< A scaling factor that the field is multiplied by before
2589		!! it is written [a A-1 ~> 1], for example to convert it
2590		!! from its internal units to the desired units for output.
2591		!! Here scale and unscale are synonymous, but unscale
2592		!! takes precedence if both are present.
2593		logical, optional, intent(in) :: zero_zeros !< If present and true, convert negative zeros
2594		!! into ordinary signless zeros.
2595
2596		! Local variables
2597	0	real, allocatable :: field_rot(:,:,:,:) ! A rotated version of field, with the same units [a] or
2598		! rescaled [A ~> a] then [a]
2599		real :: scale_fac ! A scaling factor to use before writing the array [a A-1 ~> 1]
2600		integer :: qturns ! The number of quarter turns through which to rotate field
2601
2602	0	qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4)
2603	0	scale_fac = 1.0 ; if (present(scale)) scale_fac = scale
2604	0	if (present(unscale)) scale_fac = unscale
2605
2606	0	if ((qturns == 0) .and. (scale_fac == 1.0) .and. .not.present(zero_zeros)) then
2607		call write_field(IO_handle, field_md, MOM_domain, field, tstamp=tstamp, &
2608	0	tile_count=tile_count, fill_value=fill_value)
2609		else
2610	0	call allocate_rotated_array(field, [1,1,1,1], qturns, field_rot)
2611	0	call rotate_array(field, qturns, field_rot)
2612	0	call rescale_comp_data(MOM_Domain, field_rot, scale_fac, zero_zeros)
2613		call write_field(IO_handle, field_md, MOM_domain, field_rot, tstamp=tstamp, &
2614	0	tile_count=tile_count, fill_value=fill_value)
2615	0	deallocate(field_rot)
2616		endif
2617	0	end subroutine MOM_write_field_legacy_4d
2618
2619
2620		!> Write a 3d field to an output file, potentially with rotation
2621	0	subroutine MOM_write_field_legacy_3d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, &
2622		fill_value, turns, scale, unscale, zero_zeros)
2623		type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing
2624		type(fieldtype), intent(in) :: field_md !< Field type with metadata
2625		type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition
2626		real, dimension(:,:,:), intent(inout) :: field !< Unrotated field to write in arbitrary units [A ~> a]
2627		real, optional, intent(in) :: tstamp !< Model timestamp, often in [days]
2628		integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1)
2629		real, optional, intent(in) :: fill_value !< Missing data fill value in the units used in the file [a]
2630		integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data
2631		real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied by before
2632		!! it is written [a A-1 ~> 1], for example to convert it
2633		!! from its internal units to the desired units for output
2634		real, optional, intent(in) :: unscale !< A scaling factor that the field is multiplied by before
2635		!! it is written [a A-1 ~> 1], for example to convert it
2636		!! from its internal units to the desired units for output.
2637		!! Here scale and unscale are synonymous, but unscale
2638		!! takes precedence if both are present.
2639		logical, optional, intent(in) :: zero_zeros !< If present and true, convert negative zeros
2640		!! into ordinary signless zeros.
2641
2642		! Local variables
2643	0	real, allocatable :: field_rot(:,:,:) ! A rotated version of field, with the same units [a] or
2644		! rescaled [A ~> a] then [a]
2645		real :: scale_fac ! A scaling factor to use before writing the array [a A-1 ~> 1]
2646		integer :: qturns ! The number of quarter turns through which to rotate field
2647
2648	0	qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4)
2649	0	scale_fac = 1.0 ; if (present(scale)) scale_fac = scale
2650	0	if (present(unscale)) scale_fac = unscale
2651
2652	0	if ((qturns == 0) .and. (scale_fac == 1.0) .and. .not.present(zero_zeros)) then
2653		call write_field(IO_handle, field_md, MOM_domain, field, tstamp=tstamp, &
2654	0	tile_count=tile_count, fill_value=fill_value)
2655		else
2656	0	call allocate_rotated_array(field, [1,1,1], qturns, field_rot)
2657	0	call rotate_array(field, qturns, field_rot)
2658	0	call rescale_comp_data(MOM_Domain, field_rot, scale_fac, zero_zeros)
2659		call write_field(IO_handle, field_md, MOM_domain, field_rot, tstamp=tstamp, &
2660	0	tile_count=tile_count, fill_value=fill_value)
2661	0	deallocate(field_rot)
2662		endif
2663	0	end subroutine MOM_write_field_legacy_3d
2664
2665
2666		!> Write a 2d field to an output file, potentially with rotation
2667	0	subroutine MOM_write_field_legacy_2d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, &
2668		fill_value, turns, scale, unscale, zero_zeros)
2669		type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing
2670		type(fieldtype), intent(in) :: field_md !< Field type with metadata
2671		type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition
2672		real, dimension(:,:), intent(inout) :: field !< Unrotated field to write in arbitrary units [A ~> a]
2673		real, optional, intent(in) :: tstamp !< Model timestamp, often in [days]
2674		integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1)
2675		real, optional, intent(in) :: fill_value !< Missing data fill value
2676		integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data
2677		real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied by before
2678		!! it is written [a A-1 ~> 1], for example to convert it
2679		!! from its internal units to the desired units for output
2680		real, optional, intent(in) :: unscale !< A scaling factor that the field is multiplied by before
2681		!! it is written [a A-1 ~> 1], for example to convert it
2682		!! from its internal units to the desired units for output.
2683		!! Here scale and unscale are synonymous, but unscale
2684		!! takes precedence if both are present.
2685		logical, optional, intent(in) :: zero_zeros !< If present and true, convert negative zeros
2686		!! into ordinary signless zeros.
2687
2688		! Local variables
2689	0	real, allocatable :: field_rot(:,:) ! A rotated version of field, with the same units [a] or
2690		! rescaled [A ~> a] then [a]
2691		real :: scale_fac ! A scaling factor to use before writing the array [a A-1 ~> 1]
2692		integer :: qturns ! The number of quarter turns through which to rotate field
2693
2694	0	qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4)
2695	0	scale_fac = 1.0 ; if (present(scale)) scale_fac = scale
2696	0	if (present(unscale)) scale_fac = unscale
2697
2698	0	if ((qturns == 0) .and. (scale_fac == 1.0) .and. .not.present(zero_zeros)) then
2699		call write_field(IO_handle, field_md, MOM_domain, field, tstamp=tstamp, &
2700	0	tile_count=tile_count, fill_value=fill_value)
2701		else
2702	0	call allocate_rotated_array(field, [1,1], qturns, field_rot)
2703	0	call rotate_array(field, qturns, field_rot)
2704	0	call rescale_comp_data(MOM_Domain, field_rot, scale_fac, zero_zeros)
2705		call write_field(IO_handle, field_md, MOM_domain, field_rot, tstamp=tstamp, &
2706	0	tile_count=tile_count, fill_value=fill_value)
2707	0	deallocate(field_rot)
2708		endif
2709	0	end subroutine MOM_write_field_legacy_2d
2710
2711
2712		!> Write a 1d field to an output file
2713	0	subroutine MOM_write_field_legacy_1d(IO_handle, field_md, field, tstamp, fill_value, scale, unscale, zero_zeros)
2714		type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing
2715		type(fieldtype), intent(in) :: field_md !< Field type with metadata
2716		real, dimension(:), intent(in) :: field !< Field to write in arbitrary units [A ~> a]
2717		real, optional, intent(in) :: tstamp !< Model timestamp, often in [days]
2718		real, optional, intent(in) :: fill_value !< Missing data fill value [a]
2719		real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied by before
2720		!! it is written [a A-1 ~> 1], for example to convert it
2721		!! from its internal units to the desired units for output
2722		real, optional, intent(in) :: unscale !< A scaling factor that the field is multiplied by before
2723		!! it is written [a A-1 ~> 1], for example to convert it
2724		!! from its internal units to the desired units for output.
2725		!! Here scale and unscale are synonymous, but unscale
2726		!! takes precedence if both are present.
2727		logical, optional, intent(in) :: zero_zeros !< If present and true, convert negative zeros
2728		!! into ordinary signless zeros.
2729
2730		! Local variables
2731	0	real, dimension(:), allocatable :: array ! A rescaled copy of field [a]
2732		real :: scale_fac ! A scaling factor to use before writing the array [a A-1 ~> 1]
2733		logical :: design_zeros ! If true, convert negative zeros into ordinary signless zeros.
2734		integer :: i
2735
2736	0	scale_fac = 1.0 ; if (present(scale)) scale_fac = scale
2737	0	if (present(unscale)) scale_fac = unscale
2738
2739	0	design_zeros = .false. ; if (present(zero_zeros)) design_zeros = zero_zeros
2740
2741	0	if ((scale_fac == 1.0) .and. (.not.design_zeros)) then
2742	0	call write_field(IO_handle, field_md, field, tstamp=tstamp)
2743		else
2744	0	allocate(array(size(field)))
2745	0	array(:) = scale_fac * field(:)
2746	0	if (present(fill_value)) then
2747	0	do i=1,size(field) ; if (field(i) == fill_value) array(i) = fill_value ; enddo
2748		endif
2749	0	if (design_zeros) then ! Convert negative zeros into zeros
2750	0	do i=1,size(field) ; if (array(i) == 0.0) array(i) = 0.0 ; enddo
2751		endif
2752	0	call write_field(IO_handle, field_md, array, tstamp=tstamp)
2753	0	deallocate(array)
2754		endif
2755	0	end subroutine MOM_write_field_legacy_1d
2756
2757
2758		!> Write a 0d field to an output file
2759	0	subroutine MOM_write_field_legacy_0d(IO_handle, field_md, field, tstamp, fill_value, scale, unscale, zero_zeros)
2760		type(file_type), intent(inout) :: IO_handle !< Handle for a file that is open for writing
2761		type(fieldtype), intent(in) :: field_md !< Field type with metadata
2762		real, intent(in) :: field !< Field to write in arbitrary units [A ~> a]
2763		real, optional, intent(in) :: tstamp !< Model timestamp, often in [days]
2764		real, optional, intent(in) :: fill_value !< Missing data fill value [a]
2765		real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied by before
2766		!! it is written [a A-1 ~> 1], for example to convert it
2767		!! from its internal units to the desired units for output
2768		real, optional, intent(in) :: unscale !< A scaling factor that the field is multiplied by before
2769		!! it is written [a A-1 ~> 1], for example to convert it
2770		!! from its internal units to the desired units for output.
2771		!! Here scale and unscale are synonymous, but unscale
2772		!! takes precedence if both are present.
2773		logical, optional, intent(in) :: zero_zeros !< If present and true, convert negative zeros
2774		!! into ordinary signless zeros.
2775
2776		! Local variables
2777		real :: scale_fac ! A scaling factor to use before writing the field [a A-1 ~> 1]
2778		real :: scaled_val ! A rescaled copy of field [a]
2779
2780	0	scale_fac = 1.0 ; if (present(scale)) scale_fac = scale
2781	0	if (present(unscale)) scale_fac = unscale
2782
2783	0	scaled_val = field * scale_fac
2784
2785	0	if (present(fill_value)) then ; if (field == fill_value) scaled_val = fill_value ; endif
2786	0	if (present(zero_zeros)) then ; if (zero_zeros .and. (scaled_val == 0.0)) scaled_val = 0.0 ; endif
2787
2788	0	call write_field(IO_handle, field_md, scaled_val, tstamp=tstamp)
2789	0	end subroutine MOM_write_field_legacy_0d
2790
2791
2792		!> Write a 4d field to an output file, potentially with rotation
2793	0	subroutine MOM_write_field_4d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, &
2794		fill_value, turns, scale, unscale, zero_zeros)
2795		class(MOM_file), intent(inout) :: IO_handle !< Handle for a file that is open for writing
2796		type(MOM_field), intent(in) :: field_md !< Field type with metadata
2797		type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition
2798		real, dimension(:,:,:,:), intent(inout) :: field !< Unrotated field to write in arbitrary units [A ~> a]
2799		real, optional, intent(in) :: tstamp !< Model timestamp, often in [days]
2800		integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1)
2801		real, optional, intent(in) :: fill_value !< Missing data fill value [a]
2802		integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data
2803		real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied by before
2804		!! it is written [a A-1 ~> 1], for example to convert it
2805		!! from its internal units to the desired units for output
2806		real, optional, intent(in) :: unscale !< A scaling factor that the field is multiplied by before
2807		!! it is written [a A-1 ~> 1], for example to convert it
2808		!! from its internal units to the desired units for output.
2809		!! Here scale and unscale are synonymous, but unscale
2810		!! takes precedence if both are present.
2811		logical, optional, intent(in) :: zero_zeros !< If present and true, convert negative zeros
2812		!! into ordinary signless zeros.
2813
2814		! Local variables
2815	0	real, allocatable :: field_rot(:,:,:,:) ! A rotated version of field, with the same units or rescaled [a]
2816		real :: scale_fac ! A scaling factor to use before writing the array [a A-1 ~> 1]
2817		integer :: qturns ! The number of quarter turns through which to rotate field
2818
2819	0	qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4)
2820	0	scale_fac = 1.0 ; if (present(scale)) scale_fac = scale
2821	0	if (present(unscale)) scale_fac = unscale
2822
2823	0	if ((qturns == 0) .and. (scale_fac == 1.0) .and. .not.present(zero_zeros)) then
2824		call IO_handle%write_field(field_md, MOM_domain, field, tstamp=tstamp, &
2825	0	tile_count=tile_count, fill_value=fill_value)
2826		else
2827	0	call allocate_rotated_array(field, [1,1,1,1], qturns, field_rot)
2828	0	call rotate_array(field, qturns, field_rot)
2829	0	call rescale_comp_data(MOM_Domain, field_rot, scale_fac, zero_zeros)
2830		call IO_handle%write_field(field_md, MOM_domain, field_rot, tstamp=tstamp, &
2831	0	tile_count=tile_count, fill_value=fill_value)
2832	0	deallocate(field_rot)
2833		endif
2834	0	end subroutine MOM_write_field_4d
2835
2836		!> Write a 3d field to an output file, potentially with rotation
2837	31	subroutine MOM_write_field_3d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, &
2838		fill_value, turns, scale, unscale, zero_zeros)
2839		class(MOM_file), intent(inout) :: IO_handle !< Handle for a file that is open for writing
2840		type(MOM_field), intent(in) :: field_md !< Field type with metadata
2841		type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition
2842		real, dimension(:,:,:), intent(inout) :: field !< Unrotated field to write in arbitrary units [A ~> a]
2843		real, optional, intent(in) :: tstamp !< Model timestamp, often in [days]
2844		integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1)
2845		real, optional, intent(in) :: fill_value !< Missing data fill value [a]
2846		integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data
2847		real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied by before
2848		!! it is written [a A-1 ~> 1], for example to convert it
2849		!! from its internal units to the desired units for output
2850		real, optional, intent(in) :: unscale !< A scaling factor that the field is multiplied by before
2851		!! it is written [a A-1 ~> 1], for example to convert it
2852		!! from its internal units to the desired units for output.
2853		!! Here scale and unscale are synonymous, but unscale
2854		!! takes precedence if both are present.
2855		logical, optional, intent(in) :: zero_zeros !< If present and true, convert negative zeros
2856		!! into ordinary signless zeros.
2857
2858		! Local variables
2859	31	real, allocatable :: field_rot(:,:,:) ! A rotated version of field, with the same units or rescaled [a]
2860		real :: scale_fac ! A scaling factor to use before writing the array [a A-1 ~> 1]
2861		integer :: qturns ! The number of quarter turns through which to rotate field
2862
2863	31	qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4)
2864	31	scale_fac = 1.0 ; if (present(scale)) scale_fac = scale
2865	31	if (present(unscale)) scale_fac = unscale
2866
2867	31	if ((qturns == 0) .and. (scale_fac == 1.0) .and. .not.present(zero_zeros)) then
2868		call IO_handle%write_field(field_md, MOM_domain, field, tstamp=tstamp, &
2869	0	tile_count=tile_count, fill_value=fill_value)
2870		else
2871	31	call allocate_rotated_array(field, [1,1,1], qturns, field_rot)
2872	31	call rotate_array(field, qturns, field_rot)
2873	31	call rescale_comp_data(MOM_Domain, field_rot, scale_fac, zero_zeros)
2874		call IO_handle%write_field(field_md, MOM_domain, field_rot, tstamp=tstamp, &
2875	31	tile_count=tile_count, fill_value=fill_value)
2876	31	deallocate(field_rot)
2877		endif
2878	31	end subroutine MOM_write_field_3d
2879
2880		!> Write a 2d field to an output file, potentially with rotation
2881	28	subroutine MOM_write_field_2d(IO_handle, field_md, MOM_domain, field, tstamp, tile_count, &
2882		fill_value, turns, scale, unscale, zero_zeros)
2883		class(MOM_file), intent(inout) :: IO_handle !< Handle for a file that is open for writing
2884		type(MOM_field), intent(in) :: field_md !< Field type with metadata
2885		type(MOM_domain_type), intent(in) :: MOM_domain !< The MOM_Domain that describes the decomposition
2886		real, dimension(:,:), intent(inout) :: field !< Unrotated field to write in arbitrary units [A ~> a]
2887		real, optional, intent(in) :: tstamp !< Model timestamp, often in [days]
2888		integer, optional, intent(in) :: tile_count !< PEs per tile (default: 1)
2889		real, optional, intent(in) :: fill_value !< Missing data fill value [a]
2890		integer, optional, intent(in) :: turns !< Number of quarter-turns to rotate the data
2891		real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied by before
2892		!! it is written [a A-1 ~> 1], for example to convert it
2893		!! from its internal units to the desired units for output
2894		real, optional, intent(in) :: unscale !< A scaling factor that the field is multiplied by before
2895		!! it is written [a A-1 ~> 1], for example to convert it
2896		!! from its internal units to the desired units for output.
2897		!! Here scale and unscale are synonymous, but unscale
2898		!! takes precedence if both are present.
2899		logical, optional, intent(in) :: zero_zeros !< If present and true, convert negative zeros
2900		!! into ordinary signless zeros.
2901
2902		! Local variables
2903	28	real, allocatable :: field_rot(:,:) ! A rotated version of field, with the same units or rescaled [a]
2904		real :: scale_fac ! A scaling factor to use before writing the array [a A-1 ~> 1]
2905		integer :: qturns ! The number of quarter turns through which to rotate field
2906
2907	28	qturns = 0 ; if (present(turns)) qturns = modulo(turns, 4)
2908	28	scale_fac = 1.0 ; if (present(scale)) scale_fac = scale
2909	28	if (present(unscale)) scale_fac = unscale
2910
2911	28	if ((qturns == 0) .and. (scale_fac == 1.0) .and. .not.present(zero_zeros)) then
2912		call IO_handle%write_field(field_md, MOM_domain, field, tstamp=tstamp, &
2913	0	tile_count=tile_count, fill_value=fill_value)
2914		else
2915	28	call allocate_rotated_array(field, [1,1], qturns, field_rot)
2916	28	call rotate_array(field, qturns, field_rot)
2917	28	call rescale_comp_data(MOM_Domain, field_rot, scale_fac, zero_zeros)
2918		call IO_handle%write_field(field_md, MOM_domain, field_rot, tstamp=tstamp, &
2919	28	tile_count=tile_count, fill_value=fill_value)
2920	28	deallocate(field_rot)
2921		endif
2922	28	end subroutine MOM_write_field_2d
2923
2924		!> Write a 1d field to an output file
2925	2	subroutine MOM_write_field_1d(IO_handle, field_md, field, tstamp, fill_value, scale, unscale, zero_zeros)
2926		class(MOM_file), intent(inout) :: IO_handle !< Handle for a file that is open for writing
2927		type(MOM_field), intent(in) :: field_md !< Field type with metadata
2928		real, dimension(:), intent(in) :: field !< Field to write in arbitrary units [A ~> a]
2929		real, optional, intent(in) :: tstamp !< Model timestamp, often in [days]
2930		real, optional, intent(in) :: fill_value !< Missing data fill value [a]
2931		real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied by before
2932		!! it is written [a A-1 ~> 1], for example to convert it
2933		!! from its internal units to the desired units for output
2934		real, optional, intent(in) :: unscale !< A scaling factor that the field is multiplied by before
2935		!! it is written [a A-1 ~> 1], for example to convert it
2936		!! from its internal units to the desired units for output.
2937		!! Here scale and unscale are synonymous, but unscale
2938		!! takes precedence if both are present.
2939		logical, optional, intent(in) :: zero_zeros !< If present and true, convert negative zeros
2940		!! into ordinary signless zeros.
2941
2942		! Local variables
2943	2	real, dimension(:), allocatable :: array ! A rescaled copy of field in arbtrary unscaled units [a]
2944		real :: scale_fac ! A scaling factor to use before writing the array [a A-1 ~> 1]
2945		logical :: design_zeros ! If true, convert negative zeros into ordinary signless zeros.
2946		integer :: i
2947
2948	2	scale_fac = 1.0 ; if (present(scale)) scale_fac = scale
2949	2	if (present(unscale)) scale_fac = unscale
2950
2951	2	design_zeros = .false. ; if (present(zero_zeros)) design_zeros = zero_zeros
2952
2953	2	if ((scale_fac == 1.0) .and. (.not.design_zeros)) then
2954	2	call IO_handle%write_field(field_md, field, tstamp=tstamp)
2955		else
2956	0	allocate(array(size(field)))
2957	0	array(:) = scale_fac * field(:)
2958	0	if (present(fill_value)) then
2959	0	do i=1,size(field) ; if (field(i) == fill_value) array(i) = fill_value ; enddo
2960		endif
2961	0	if (design_zeros) then ! Convert negative zeros into zeros
2962	0	do i=1,size(field) ; if (array(i) == 0.0) array(i) = 0.0 ; enddo
2963		endif
2964	0	call IO_handle%write_field(field_md, array, tstamp=tstamp)
2965	0	deallocate(array)
2966		endif
2967	2	end subroutine MOM_write_field_1d
2968
2969		!> Write a 0d field to an output file
2970	4	subroutine MOM_write_field_0d(IO_handle, field_md, field, tstamp, fill_value, scale, unscale, zero_zeros)
2971		class(MOM_file), intent(inout) :: IO_handle !< Handle for a file that is open for writing
2972		type(MOM_field), intent(in) :: field_md !< Field type with metadata
2973		real, intent(in) :: field !< Field to write in arbitrary units [A ~> a]
2974		real, optional, intent(in) :: tstamp !< Model timestamp, often in [days]
2975		real, optional, intent(in) :: fill_value !< Missing data fill value [a]
2976		real, optional, intent(in) :: scale !< A scaling factor that the field is multiplied by before
2977		!! it is written [a A-1 ~> 1], for example to convert it
2978		!! from its internal units to the desired units for output
2979		real, optional, intent(in) :: unscale !< A scaling factor that the field is multiplied by before
2980		!! it is written [a A-1 ~> 1], for example to convert it
2981		!! from its internal units to the desired units for output.
2982		!! Here scale and unscale are synonymous, but unscale
2983		!! takes precedence if both are present.
2984		logical, optional, intent(in) :: zero_zeros !< If present and true, convert negative zeros
2985		!! into ordinary signless zeros.
2986
2987		! Local variables
2988		real :: scale_fac ! A scaling factor to use before writing the field [a A-1 ~> 1]
2989		real :: scaled_val ! A rescaled copy of field in arbtrary unscaled units [a]
2990
2991	4	scale_fac = 1.0 ; if (present(scale)) scale_fac = scale
2992	4	if (present(unscale)) scale_fac = unscale
2993
2994	4	scaled_val = field * scale_fac
2995
2996	4	if (present(fill_value)) then ; if (field == fill_value) scaled_val = fill_value ; endif
2997	4	if (present(zero_zeros)) then ; if (zero_zeros .and. (scaled_val == 0.0)) scaled_val = 0.0 ; endif
2998
2999	4	call IO_handle%write_field(field_md, scaled_val, tstamp=tstamp)
3000	4	end subroutine MOM_write_field_0d
3001
3002		!> Given filename and fieldname, this subroutine returns the size of the field in the file
3003	0	subroutine field_size(filename, fieldname, sizes, field_found, no_domain, ndims, ncid_in)
3004		character(len=*), intent(in) :: filename !< The name of the file to read
3005		character(len=*), intent(in) :: fieldname !< The name of the variable whose sizes are returned
3006		integer, dimension(:), intent(inout) :: sizes !< The sizes of the variable in each dimension
3007		logical, optional, intent(out) :: field_found !< This indicates whether the field was found in
3008		!! the input file. Without this argument, there
3009		!! is a fatal error if the field is not found.
3010		logical, optional, intent(in) :: no_domain !< If present and true, do not check for file
3011		!! names with an appended tile number. If
3012		!! ndims is present, the default changes to true.
3013		integer, optional, intent(out) :: ndims !< The number of dimensions to the variable
3014		integer, optional, intent(in) :: ncid_in !< The netCDF ID of an open file. If absent, the
3015		!! file is opened and closed within this routine.
3016
3017	0	if (present(ndims)) then
3018	0	if (present(no_domain)) then ; if (.not.no_domain) call MOM_error(FATAL, &
3019	0	"field_size does not support the ndims argument when no_domain is present and false.")
3020		endif
3021	0	call get_var_sizes(filename, fieldname, ndims, sizes, match_case=.false., ncid_in=ncid_in)
3022	0	if (present(field_found)) field_found = (ndims >= 0)
3023	0	if ((ndims < 0) .and. .not.present(field_found)) then
3024	0	call MOM_error(FATAL, "Variable "//trim(fieldname)//" not found in "//trim(filename) )
3025		endif
3026		else
3027	0	call get_field_size(filename, fieldname, sizes, field_found=field_found, no_domain=no_domain)
3028		endif
3029
3030	0	end subroutine field_size
3031
3032
3033		!> Copies a string
3034	994	subroutine safe_string_copy(str1, str2, fieldnm, caller)
3035		character(len=*), intent(in) :: str1 !< The string being copied
3036		character(len=*), intent(out) :: str2 !< The string being copied into
3037		character(len=*), optional, intent(in) :: fieldnm !< The name of the field for error messages
3038		character(len=*), optional, intent(in) :: caller !< The calling routine for error messages
3039
3040	994	if (len(trim(str1)) > len(str2)) then
3041	0	if (present(fieldnm) .and. present(caller)) then
3042		call MOM_error(FATAL, trim(caller)//" attempted to copy the overly long string "//&
3043	0	trim(str1)//" into "//trim(fieldnm))
3044		else
3045		call MOM_error(FATAL, "safe_string_copy: The string "//trim(str1)//&
3046	0	" is longer than its intended target.")
3047		endif
3048		endif
3049	994	str2 = trim(str1)
3050	994	end subroutine safe_string_copy
3051
3052		!> Returns a name with "%#E" or "%E" replaced with the ensemble member number.
3053	11	function ensembler(name, ens_no_in) result(en_nm)
3054		character(len=*), intent(in) :: name !< The name to be modified
3055		integer, optional, intent(in) :: ens_no_in !< The number of the current ensemble member
3056		character(len=len(name)) :: en_nm !< The name encoded with the ensemble number
3057
3058		! This function replaces "%#E" or "%E" with the ensemble number anywhere it
3059		! occurs in name, with %E using 4 or 6 digits (depending on the ensemble size)
3060		! and %#E using # digits, where # is a number from 1 to 9.
3061
3062	11	character(len=len(name)) :: tmp
3063		character(10) :: ens_num_char
3064		character(3) :: code_str
3065		integer :: ens_no
3066		integer :: n, is
3067
3068	11	en_nm = trim(name)
3069	11	if (index(name,"%") == 0) return
3070
3071	0	if (present(ens_no_in)) then
3072	0	ens_no = ens_no_in
3073		else
3074	0	ens_no = get_ensemble_id()
3075		endif
3076
3077	0	write(ens_num_char, '(I0)') ens_no
3078	0	do
3079	0	is = index(en_nm,"%E")
3080	0	if (is == 0) exit
3081	0	if (len(en_nm) < len(trim(en_nm)) + len(trim(ens_num_char)) - 2) &
3082		call MOM_error(FATAL, "MOM_io ensembler: name "//trim(name)// &
3083	0	" is not long enough for %E expansion for ens_no "//trim(ens_num_char))
3084	0	tmp = en_nm(1:is-1)//trim(ens_num_char)//trim(en_nm(is+2:))
3085	0	en_nm = tmp
3086		enddo
3087
3088	0	if (index(name,"%") == 0) return
3089
3090	0	write(ens_num_char, '(I10.10)') ens_no
3091	0	do n=1,9 ; do
3092	0	write(code_str, '("%",I1,"E")') n
3093
3094	0	is = index(en_nm,code_str)
3095	0	if (is == 0) exit
3096	0	if (ens_no < 10**n) then
3097	0	if (len(en_nm) < len(trim(en_nm)) + n-3) call MOM_error(FATAL, &
3098	0	"MOM_io ensembler: name "//trim(name)//" is not long enough for %E expansion.")
3099	0	tmp = en_nm(1:is-1)//trim(ens_num_char(11-n:10))//trim(en_nm(is+3:))
3100		else
3101		call MOM_error(FATAL, "MOM_io ensembler: Ensemble number is too large "//&
3102	0	"to be encoded with "//code_str//" in "//trim(name))
3103		endif
3104	0	en_nm = tmp
3105		enddo ; enddo
3106
3107	11	end function ensembler
3108
3109		!> Provide a string to append to filenames, to differentiate ensemble members, for example.
3110	7	subroutine get_filename_appendix(suffix)
3111		character(len=*), intent(out) :: suffix !< A string to append to filenames
3112
3113	7	call get_filename_suffix(suffix)
3114	7	end subroutine get_filename_appendix
3115
3116		!> Write a file version number to the log file or other output file
3117	0	subroutine write_version_number(version, tag, unit)
3118		character(len=*), intent(in) :: version !< A string that contains the routine name and version
3119		character(len=*), optional, intent(in) :: tag !< A tag name to add to the message
3120		integer, optional, intent(in) :: unit !< An alternate unit number for output
3121
3122	0	call write_version(version, tag, unit)
3123	0	end subroutine write_version_number
3124
3125
3126		!> Open a single namelist file that is potentially readable by all PEs.
3127	2	function open_namelist_file(file) result(unit)
3128		character(len=*), optional, intent(in) :: file !< The file to open, by default "input.nml"
3129		integer :: unit !< The opened unit number of the namelist file
3130	2	unit = MOM_namelist_file(file)
3131	4	end function open_namelist_file
3132
3133		!> Checks the iostat argument that is returned after reading a namelist variable and writes a
3134		!! message if there is an error.
3135	3	function check_nml_error(IOstat, nml_name) result(ierr)
3136		integer, intent(in) :: IOstat !< An I/O status field from a namelist read call
3137		character(len=*), intent(in) :: nml_name !< The name of the namelist
3138		integer :: ierr !< A copy of IOstat that is returned to preserve legacy function behavior
3139	3	call check_namelist_error(IOstat, nml_name)
3140	3	ierr = IOstat
3141	3	end function check_nml_error
3142
3143		!> Initialize the MOM_io module
3144	1	subroutine MOM_io_init(param_file)
3145		type(param_file_type), intent(in) :: param_file !< structure indicating the open file to
3146		!! parse for model parameter values.
3147
3148		! This include declares and sets the variable "version".
3149		# include "version_variable.h"
3150		character(len=40) :: mdl = "MOM_io" ! This module's name.
3151
3152	1	call log_version(param_file, mdl, version)
3153
3154	1	end subroutine MOM_io_init
3155		!> Returns the dimension variable information for a netCDF variable
3156	0	subroutine get_var_axes_info(filename, fieldname, axes_info)
3157		character(len=*), intent(in) :: filename !< A filename from which to read
3158		character(len=*), intent(in) :: fieldname !< The name of the field to read
3159		type(axis_info), dimension(4), intent(inout) :: axes_info !< A returned array of field axis information
3160
3161		!! local variables
3162		integer :: rcode
3163		logical :: success
3164		integer :: ncid, varid, ndims
3165		integer :: id, jd, kd
3166		integer, dimension(4) :: dims, dim_id
3167		character(len=128) :: dim_name(4)
3168		integer, dimension(1) :: start, count
3169		!! cartesian axis data
3170	0	real, allocatable, dimension(:) :: x ! x-axis labels, often [degrees_E] or [km] or [m]
3171	0	real, allocatable, dimension(:) :: y ! y-axis labels, often [degrees_N] or [km] or [m]
3172	0	real, allocatable, dimension(:) :: z ! vertical axis labels [various], often [m] or [kg m-3]
3173
3174
3175	0	call open_file_to_read(filename, ncid, success=success)
3176
3177	0	rcode = NF90_INQ_VARID(ncid, trim(fieldname), varid)
3178	0	if (rcode /= 0) call MOM_error(FATAL,"error finding variable "//trim(fieldname)//&
3179	0	" in file "//trim(filename)//" in hinterp_extrap")
3180
3181	0	rcode = NF90_INQUIRE_VARIABLE(ncid, varid, ndims=ndims, dimids=dims)
3182	0	if (rcode /= 0) call MOM_error(FATAL, "Error inquiring about the dimensions of "//trim(fieldname)//&
3183	0	" in file "//trim(filename)//" in hinterp_extrap")
3184	0	if (ndims < 3) call MOM_error(FATAL,"Variable "//trim(fieldname)//" in file "//trim(filename)// &
3185	0	" has too few dimensions to be read as a 3-d array.")
3186	0	rcode = NF90_INQUIRE_DIMENSION(ncid, dims(1), dim_name(1), len=id)
3187	0	if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 1 data for "// &
3188	0	trim(fieldname)//" in file "// trim(filename)//" in hinterp_extrap")
3189	0	rcode = NF90_INQ_VARID(ncid, dim_name(1), dim_id(1))
3190	0	if (rcode /= 0) call MOM_error(FATAL,"error finding variable "//trim(dim_name(1))//&
3191	0	" in file "//trim(filename)//" in hinterp_extrap")
3192	0	rcode = NF90_INQUIRE_DIMENSION(ncid, dims(2), dim_name(2), len=jd)
3193	0	if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 2 data for "// &
3194	0	trim(fieldname)//" in file "// trim(filename)//" in hinterp_extrap")
3195	0	rcode = NF90_INQ_VARID(ncid, dim_name(2), dim_id(2))
3196	0	if (rcode /= 0) call MOM_error(FATAL,"error finding variable "//trim(dim_name(2))//&
3197	0	" in file "//trim(filename)//" in hinterp_extrap")
3198	0	rcode = NF90_INQUIRE_DIMENSION(ncid, dims(3), dim_name(3), len=kd)
3199	0	if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 3 data for "// &
3200	0	trim(fieldname)//" in file "// trim(filename)//" in hinterp_extrap")
3201	0	rcode = NF90_INQ_VARID(ncid, dim_name(3), dim_id(3))
3202	0	if (rcode /= 0) call MOM_error(FATAL,"error finding variable "//trim(dim_name(3))//&
3203	0	" in file "//trim(filename)//" in hinterp_extrap")
3204	0	allocate(x(id), y(jd), z(kd))
3205
3206	0	start = 1 ; count = 1 ; count(1) = id
3207	0	rcode = NF90_GET_VAR(ncid, dim_id(1), x, start, count)
3208	0	if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 1 values for var_name "// &
3209	0	trim(fieldname)//",dim_name "//trim(dim_name(1))//" in file "// trim(filename)//" in hinterp_extrap")
3210	0	start = 1 ; count = 1 ; count(1) = jd
3211	0	rcode = NF90_GET_VAR(ncid, dim_id(2), y, start, count)
3212	0	if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 2 values for var_name "// &
3213	0	trim(fieldname)//",dim_name "//trim(dim_name(2))//" in file "// trim(filename)//" in hinterp_extrap")
3214	0	start = 1 ; count = 1 ; count(1) = kd
3215	0	rcode = NF90_GET_VAR(ncid, dim_id(3), z, start, count)
3216	0	if (rcode /= 0) call MOM_error(FATAL,"error reading dimension 3 values for var_name "// &
3217	0	trim(fieldname//",dim_name "//trim(dim_name(3)))//" in file "// trim(filename)//" in hinterp_extrap")
3218
3219	0	call set_axis_info(axes_info(1), name=trim(dim_name(1)), ax_size=id, ax_data=x,cartesian='X')
3220	0	call set_axis_info(axes_info(2), name=trim(dim_name(2)), ax_size=jd, ax_data=y,cartesian='Y')
3221	0	call set_axis_info(axes_info(3), name=trim(dim_name(3)), ax_size=kd, ax_data=z,cartesian='Z')
3222
3223	0	call close_file_to_read(ncid, filename)
3224
3225	0	deallocate(x,y,z)
3226
3227	0	end subroutine get_var_axes_info
3228		!> \namespace mom_io
3229		!!
3230		!! This file contains a number of subroutines that manipulate
3231		!! NetCDF files and handle input and output of fields. These
3232		!! subroutines, along with their purpose, are:
3233		!!
3234		!! * create_file: create a new file and set up structures that are
3235		!! needed for subsequent output and write out the coordinates.
3236		!! * reopen_file: reopen an existing file for writing and set up
3237		!! structures that are needed for subsequent output.
3238		!! * open_input_file: open the indicated file for reading only.
3239		!! * close_file: close an open file.
3240		!! * synch_file: flush the buffers, completing all pending output.
3241		!!
3242		!! * write_field: write a field to an open file.
3243		!! * write_time: write a value of the time axis to an open file.
3244		!! * read_data: read a variable from an open file.
3245		!! * read_time: read a time from an open file.
3246		!!
3247		!! * name_output_file: provide a name for an output file based on a
3248		!! name root and the time of the output.
3249		!! * find_input_file: find a file that has been previously written by
3250		!! MOM and named by name_output_file and open it for reading.
3251		!!
3252		!! * handle_error: write an error code and quit.
3253
3254	0	end module MOM_io