Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

mpp_io_mod.F90 @ 1677

Last change on this file since 1677 was 1677, checked in by aclsce, 12 years ago
Imported oasis3 (tag ipslcm5a) from cvs server to svn server (igcmg project).
File size: 137.9 KB

Line
1	!-----------------------------------------------------------------------
2	! Parallel I/O for message-passing codes
3	!
4	! AUTHOR: V. Balaji (vbalaji@noaa.gov)
5	! Princeton University/GFDL
6	!
7	! MODIFICATIONS: Reiner Vogelsang (reiner@sgi.com)
8	!
9	! This program is free software; The author agrees that you can
10	! redistribute and/or modify this version of the program under the
11	! terms of the Lesser GNU General Public License as published
12	! by the Free Software Foundation.
13	!
14	! This program is distributed in the hope that it will be useful,
15	! but WITHOUT ANY WARRANTY; without even the implied warranty of
16	! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17	! Lesser GNU General Public License for more details
18	! (http://www.gnu.org/copyleft/lesser.html).
19	!-----------------------------------------------------------------------
20	!#include <os.h>
21
22	module mpp_io_mod
23	use mod_kinds_model
24	use mpp_mod
25	use mpp_domains_mod
26	implicit none
27	#include <os.h>
28	private
29
30	character(len=128), private :: version= &
31	'$Id: mpp_io_mod.F90,v 1.1.1.1 2005/03/23 16:01:11 adm Exp $'
32	character(len=128), private :: tagname= &
33	'$Name: ipslcm5a $'
34
35	integer, private :: pe, npes
36
37	type, public :: axistype
38	private
39	character(len=128) :: name
40	character(len=128) :: units
41	character(len=256) :: longname
42	character(len=8) :: cartesian
43	integer :: sense, len !+/-1, depth or height?
44	type(domain1D) :: domain !if pointer is associated, it is a distributed data axis
45	real, pointer :: data(:) !axis values (not used if time axis)
46	integer :: id, did, type, natt !id is the "variable ID", did is the "dimension ID": netCDF requires 2 IDs for axes
47	type(atttype), pointer :: Att(:)
48	end type axistype
49
50	type, public :: atttype
51	integer :: type, len
52	character(len=128) :: name
53	character(len=256) :: catt
54	! just use type conversion for integers
55	real, pointer :: fatt(:)
56	end type atttype
57
58	type, public :: fieldtype
59	private
60	character(len=128) :: name
61	character(len=128) :: units
62	character(len=256) :: longname
63	real :: min, max, missing, fill, scale, add
64	integer :: pack
65	type(axistype), pointer :: axes(:) !axes associated with field
66	!size, time_axis_index redundantly hold info already contained in axes
67	!it's clunky and inelegant, but required so that axes can be shared among multiple files
68	integer, pointer :: size(:)
69	integer :: time_axis_index
70	integer :: id, type, natt, ndim
71	type(atttype), pointer :: Att(:)
72	end type fieldtype
73
74	type, private :: filetype
75	character(len=256) :: name
76	integer :: action, format, access, threading, fileset, record, ncid
77	logical :: opened, initialized, nohdrs
78	integer :: time_level
79	real(DOUBLE_KIND) :: time
80	integer :: id !variable ID of time axis associated with file (only one time axis per file)
81	integer :: recdimid !dim ID of time axis associated with file (only one time axis per file)
82	!
83	! time axis values are stored here instead of axis%data since mpp_write
84	! assumes these values are not time values. Not used in mpp_write
85	!
86	real(DOUBLE_KIND), pointer :: time_values(:)
87
88	! additional elements of filetype for mpp_read (ignored for mpp_write)
89	integer :: ndim, nvar, natt ! number of dimensions, non-dimension variables and global attributes
90	! redundant axis types stored here and in associated fieldtype
91	! some axes are not used by any fields, i.e. "edges"
92	type(axistype), pointer :: axis(:)
93	type(fieldtype), pointer :: var(:)
94	type(atttype), pointer :: att(:)
95	end type filetype
96
97	type(axistype), public :: default_axis !provided to users with default components
98	type(fieldtype), public :: default_field !provided to users with default components
99	type(atttype), public :: default_att !provided to users with default components
100	!action on open
101	integer, parameter, public :: MPP_WRONLY=100, MPP_RDONLY=101, MPP_APPEND=102, MPP_OVERWR=103
102	!format
103	integer, parameter, public :: MPP_ASCII=200, MPP_IEEE32=201, MPP_NATIVE=202, MPP_NETCDF=203
104	!access
105	integer, parameter, public :: MPP_SEQUENTIAL=300, MPP_DIRECT=301
106	!threading, fileset
107	integer, parameter, public :: MPP_SINGLE=400, MPP_MULTI=401
108	!action on close
109	integer, parameter, public :: MPP_DELETE=501, MPP_COLLECT=502
110
111	type(filetype), private, allocatable :: mpp_file(:)
112	integer, private :: records_per_pe
113	integer, private :: maxunits, unit_begin, unit_end
114	integer, private :: varnum=0
115	integer, private :: error
116	character(len=256) :: text
117	!null unit: returned by PEs not participating in IO after a collective call
118	integer, parameter, private :: NULLUNIT=-1
119	real(DOUBLE_KIND), parameter, private :: NULLTIME=-1.
120	logical, private :: verbose=.FALSE., debug=.FALSE., module_is_initialized=.FALSE.
121
122	real(DOUBLE_KIND), private, allocatable :: mpp_io_stack(:)
123	integer, private :: mpp_io_stack_size=0, mpp_io_stack_hwm=0
124
125	interface mpp_write_meta
126	module procedure mpp_write_meta_var
127	module procedure mpp_write_meta_scalar_r
128	module procedure mpp_write_meta_scalar_i
129	module procedure mpp_write_meta_axis
130	module procedure mpp_write_meta_field
131	module procedure mpp_write_meta_global
132	module procedure mpp_write_meta_global_scalar_r
133	module procedure mpp_write_meta_global_scalar_i
134	end interface
135
136	interface mpp_copy_meta
137	module procedure mpp_copy_meta_axis
138	module procedure mpp_copy_meta_field
139	module procedure mpp_copy_meta_global
140	end interface
141
142	interface mpp_write
143	module procedure mpp_write_2ddecomp_r1d
144	module procedure mpp_write_2ddecomp_r2d
145	module procedure mpp_write_2ddecomp_r3d
146	module procedure mpp_write_r0D
147	module procedure mpp_write_r1D
148	module procedure mpp_write_r2D
149	module procedure mpp_write_r3D
150	module procedure mpp_write_r4D
151	module procedure mpp_write_axis
152	end interface
153
154	interface mpp_read
155	module procedure mpp_read_2ddecomp_r1d
156	module procedure mpp_read_2ddecomp_r2d
157	module procedure mpp_read_2ddecomp_r3d
158	module procedure mpp_read_r0D
159	module procedure mpp_read_r1D
160	module procedure mpp_read_r2D
161	module procedure mpp_read_r3D
162	end interface
163
164	interface mpp_get_id
165	module procedure mpp_get_axis_id
166	module procedure mpp_get_field_id
167	end interface
168
169	interface mpp_get_atts
170	module procedure mpp_get_global_atts
171	module procedure mpp_get_field_atts
172	module procedure mpp_get_axis_atts
173	end interface
174
175	interface mpp_modify_meta
176	! module procedure mpp_modify_att_meta
177	module procedure mpp_modify_field_meta
178	module procedure mpp_modify_axis_meta
179	end interface
180
181	public :: mpp_close, mpp_flush, mpp_get_iospec, mpp_get_id, mpp_get_ncid, mpp_get_unit_range, mpp_io_init, mpp_io_exit, &
182	mpp_open, mpp_set_unit_range, mpp_write, mpp_write_meta, mpp_read, mpp_get_info, mpp_get_atts, &
183	mpp_get_fields, mpp_get_times, mpp_get_axes, mpp_copy_meta, mpp_get_recdimid, mpp_get_axis_data, mpp_modify_meta, &
184	mpp_io_set_stack_size, mpp_get_field_index
185
186	private :: read_record, mpp_read_meta, lowercase
187
188	#ifdef use_netCDF
189	#include <netcdf.inc>
190	#endif
191
192	contains
193
194	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
195	! !
196	! mpp_io_init: initialize parallel I/O !
197	! !
198	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
199	subroutine mpp_io_init( flags, maxunit,maxresunit )
200	integer, intent(in), optional :: flags, maxunit ,maxresunit
201	!rv
202	!rv I introduced the variable to indentify that the top max_reserved_units
203	!rv of maxunits are reserved for OASIS coupler specific things like the trace
204	!rv files. This variable is active only if one specifies explicitely the
205	!rv argument maxunit.
206	integer::max_reserved_units
207	!rv
208	!initialize IO package: initialize mpp_file array, set valid range of units for fortran IO
209
210	if( module_is_initialized )return
211	call mpp_init(flags) !if mpp_init has been called, this call will merely return
212	pe = mpp_pe()
213	npes = mpp_npes()
214	call mpp_domains_init(flags)
215
216	maxunits = 64
217	if( PRESENT(maxunit) )maxunits = maxunit
218
219	max_reserved_units=5
220	if( PRESENT(maxresunit) )max_reserved_units = maxresunit
221
222	if( PRESENT(flags) )then
223	debug = flags.EQ.MPP_DEBUG
224	verbose = flags.EQ.MPP_VERBOSE .OR. debug
225	end if
226	!initialize default_field
227	default_field%name = 'noname'
228	default_field%units = 'nounits'
229	default_field%longname = 'noname'
230	default_field%id = -1
231	default_field%type = -1
232	default_field%natt = -1
233	default_field%ndim = -1
234	!largest possible 4-byte reals
235	default_field%min = -huge(1._4)
236	default_field%max = huge(1._4)
237	default_field%missing = -1e36
238	default_field%fill = -1e36
239	default_field%scale = 0.
240	default_field%add = huge(1._4)
241	default_field%pack = 1
242	default_field%time_axis_index = -1 !this value will never match any index
243	! Initialize default axis
244	default_axis%name = 'noname'
245	default_axis%units = 'nounits'
246	default_axis%longname = 'noname'
247	default_axis%cartesian = 'none'
248	default_axis%sense = 0
249	default_axis%len = -1
250	default_axis%id = -1
251	default_axis%did = -1
252	default_axis%type = -1
253	default_axis%natt = -1
254	! Initialize default attribute
255	default_att%name = 'noname'
256	default_att%type = -1
257	default_att%len = -1
258	default_att%catt = 'none'
259
260	!up to MAXUNITS fortran units and MAXUNITS netCDF units are supported
261	!file attributes (opened, format, access, threading, fileset) are saved against the unit number
262	!external handles to netCDF units are saved from maxunits+1:2*maxunits
263	allocate( mpp_file(NULLUNIT:2*maxunits) ) !starts at NULLUNIT=-1, used by non-participant PEs in single-threaded I/O
264	mpp_file(:)%name = ' '
265	mpp_file(:)%action = -1
266	mpp_file(:)%format = -1
267	mpp_file(:)%threading = -1
268	mpp_file(:)%fileset = -1
269	mpp_file(:)%record = -1
270	mpp_file(:)%ncid = -1
271	mpp_file(:)%opened = .FALSE.
272	mpp_file(:)%initialized = .FALSE.
273	mpp_file(:)%time_level = 0
274	mpp_file(:)%time = NULLTIME
275	mpp_file(:)%id = -1
276	!
277	mpp_file(:)%ndim = -1
278	mpp_file(:)%nvar = -1
279	!NULLUNIT "file" is always single-threaded, open and initialized (to pass checks in mpp_write)
280	mpp_file(NULLUNIT)%threading = MPP_SINGLE
281	mpp_file(NULLUNIT)%opened = .TRUE.
282	mpp_file(NULLUNIT)%initialized = .TRUE.
283	!declare the stdunits to be open
284	mpp_file(stdin ())%opened = .TRUE.
285	mpp_file(stdout())%opened = .TRUE.
286	mpp_file(stderr())%opened = .TRUE.
287	mpp_file(stdlog())%opened = .TRUE.
288	!set range of allowed fortran unit numbers: could be compiler-dependent (should not overlap stdin/out/err)
289	!
290	!rv For OASIS 3 I consider the top max_reserved_units to be excluded from
291	!rv the list of files ito closed during mpp_io_exit.
292	!rv call mpp_set_unit_range( 7, maxunits )
293	if(present(maxunit)) then
294	call mpp_set_unit_range( 7, maxunits-max_reserved_units )
295	else
296	call mpp_set_unit_range( 7, maxunits )
297	endif
298	!rv
299
300	if( pe.EQ.mpp_root_pe() )then
301	write( stdlog(),'(/a)' )'MPP_IO module '//trim(version)
302	#ifdef use_netCDF
303	text = NF_INQ_LIBVERS()
304	write( stdlog(),'(a)' )'Using netCDF library version '//trim(text)
305	#endif
306	endif
307
308	#ifdef CRAYPVP
309	!we require every file to be assigned threadwise: PVPs default to global, and are reset here
310	call ASSIGN( 'assign -P thread p:%', error )
311	#endif
312
313	call mpp_io_set_stack_size(131072) ! default initial value
314	call mpp_sync()
315	module_is_initialized = .TRUE.
316	return
317	end subroutine mpp_io_init
318
319	subroutine mpp_io_exit()
320	integer :: unit
321
322	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_IO_EXIT: must first call mpp_io_init.' )
323	!close all open fortran units
324	do unit = unit_begin,unit_end
325	if( mpp_file(unit)%opened )call FLUSH(unit)
326	end do
327	call mpp_sync()
328	do unit = unit_begin,unit_end
329	if( mpp_file(unit)%opened )close(unit)
330	end do
331	#ifdef use_netCDF
332	!close all open netCDF units
333	do unit = maxunits+1,2*maxunits
334	if( mpp_file(unit)%opened )error = NF_CLOSE(mpp_file(unit)%ncid)
335	end do
336	#endif
337
338	call mpp_max(mpp_io_stack_hwm)
339
340	if( pe.EQ.mpp_root_pe() )then
341	! write( stdout,'(/a)' )'Exiting MPP_IO module...'
342	! write( stdout,* )'MPP_IO_STACK high water mark=', mpp_io_stack_hwm
343	end if
344	deallocate(mpp_file)
345	module_is_initialized = .FALSE.
346	return
347	end subroutine mpp_io_exit
348
349	subroutine mpp_io_set_stack_size(n)
350	!set the mpp_io_stack variable to be at least n LONG words long
351	integer, intent(in) :: n
352	character(len=8) :: text
353
354	if( n.GT.mpp_io_stack_size .AND. allocated(mpp_io_stack) )deallocate(mpp_io_stack)
355	if( .NOT.allocated(mpp_io_stack) )then
356	allocate( mpp_io_stack(n) )
357	mpp_io_stack_size = n
358	write( text,'(i8)' )n
359	if( pe.EQ.mpp_root_pe() )call mpp_error( NOTE, 'MPP_IO_SET_STACK_SIZE: stack size set to '//text//'.' )
360	end if
361
362	return
363	end subroutine mpp_io_set_stack_size
364
365	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
366	! !
367	! OPENING AND CLOSING FILES: mpp_open() and mpp_close() !
368	! !
369	! mpp_open( unit, file, action, form, access, threading, & !
370	! fileset, iospec, nohdrs, recl, pelist ) !
371	! integer, intent(out) :: unit !
372	! character(len=*), intent(in) :: file !
373	! integer, intent(in), optional :: action, form, access, threading, !
374	! fileset, recl !
375	! character(len=*), intent(in), optional :: iospec !
376	! logical, intent(in), optional :: nohdrs !
377	! integer, optional, intent(in) :: pelist(:) !default ALL !
378	! !
379	! unit is intent(OUT): always _returned_by_ mpp_open() !
380	! file is the filename: REQUIRED !
381	! we append .nc to filename if it is a netCDF file !
382	! we append .<pppp> to filename if fileset is private (pppp is PE number) !
383	! iospec is a system hint for I/O organization !
384	! e.g assign(1) on SGI/Cray systems. !
385	! if nohdrs is .TRUE. headers are not written on non-netCDF writes. !
386	! nohdrs has no effect when action=MPP_RDONLY\|MPP_APPEND !
387	! or when form=MPP_NETCDF !
388	! FLAGS: !
389	! action is one of MPP_RDONLY, MPP_APPEND or MPP_WRONLY !
390	! form is one of MPP_ASCII: formatted read/write !
391	! MPP_NATIVE: unformatted read/write, no conversion !
392	! MPP_IEEE32: unformatted read/write, conversion to IEEE32 !
393	! MPP_NETCDF: unformatted read/write, conversion to netCDF !
394	! access is one of MPP_SEQUENTIAL or MPP_DIRECT (ignored for netCDF) !
395	! RECL argument is REQUIRED for direct access IO !
396	! threading is one of MPP_SINGLE or MPP_MULTI !
397	! single-threaded IO in a multi-PE run is done by PE0 !
398	! fileset is one of MPP_MULTI and MPP_SINGLE !
399	! fileset is only used for multi-threaded I/O !
400	! if all I/O PEs in <pelist> use a single fileset, !
401	! they write to the same file !
402	! if all I/O PEs in <pelist> use a multi fileset, !
403	! they each write an independent file !
404	! recl is the record length in bytes !
405	! pelist is the list of I/O PEs (currently ALL) !
406	! !
407	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
408	subroutine mpp_open( unit, file, action, form, access, threading, &
409	fileset, iospec, nohdrs, recl, pelist )
410	integer, intent(out) :: unit
411	character(len=*), intent(in) :: file
412	integer, intent(in), optional :: action, form, access, threading, &
413	fileset, recl
414	character(len=*), intent(in), optional :: iospec
415	logical, intent(in), optional :: nohdrs
416	integer, intent(in), optional :: pelist(:) !default ALL
417
418	character(len=16) :: act, acc, for, pos
419	integer :: action_flag, form_flag, access_flag, threading_flag, fileset_flag, length
420	logical :: exists
421	character(len=64) :: filespec
422	type(axistype) :: unlim !used by netCDF with mpp_append
423
424	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_OPEN: must first call mpp_io_init.' )
425	!set flags
426	action_flag = MPP_WRONLY !default
427	if( PRESENT(action) )action_flag = action
428	form_flag = MPP_ASCII
429	if( PRESENT(form) )form_flag = form
430	#ifndef use_netCDF
431	if( form_flag.EQ.MPP_NETCDF ) &
432	call mpp_error( FATAL, 'MPP_OPEN: To open a file with form=MPP_NETCDF, you must compile mpp_io with -Duse_netCDF.' )
433	#endif
434	access_flag = MPP_SEQUENTIAL
435	if( PRESENT(access) )access_flag = access
436	threading_flag = MPP_SINGLE
437	if( npes.GT.1 .AND. PRESENT(threading) )threading_flag = threading
438	fileset_flag = MPP_MULTI
439	if( PRESENT(fileset) )fileset_flag = fileset
440	if( threading_flag.EQ.MPP_SINGLE )fileset_flag = MPP_SINGLE
441
442	!get a unit number
443	if( threading_flag.EQ.MPP_SINGLE )then
444	if( pe.NE.mpp_root_pe() .AND. action_flag.NE.MPP_RDONLY )then
445	unit = NULLUNIT !PEs not participating in IO from this mpp_open() will return this value for unit
446	return
447	end if
448	end if
449	if( form_flag.EQ.MPP_NETCDF )then
450	do unit = maxunits+1,2*maxunits
451	if( .NOT.mpp_file(unit)%opened )exit
452	end do
453	if( unit.GT.2*maxunits )call mpp_error( FATAL, 'MPP_OPEN: too many open netCDF files.' )
454	else
455	do unit = unit_begin, unit_end
456	inquire( unit,OPENED=mpp_file(unit)%opened )
457	if( .NOT.mpp_file(unit)%opened )exit
458	end do
459	if( unit.GT.unit_end )call mpp_error( FATAL, 'MPP_OPEN: no available units.' )
460	end if
461
462	!get a filename
463	text = file
464	length = len(file)
465
466	if( form_flag.EQ.MPP_NETCDF.AND. file(length-2:length) /= '.nc' ) &
467	text = trim(file)//'.nc'
468
469	if( fileset_flag.EQ.MPP_MULTI )write( text,'(a,i4.4)' )trim(text)//'.', pe
470	mpp_file(unit)%name = text
471	if( verbose )print '(a,2i3,x,a,5i5)', 'MPP_OPEN: PE, unit, filename, action, format, access, threading, fileset=', &
472	pe, unit, trim(mpp_file(unit)%name), action_flag, form_flag, access_flag, threading_flag, fileset_flag
473
474	!action: read, write, overwrite, append: act and pos are ignored by netCDF
475	if( action_flag.EQ.MPP_RDONLY )then
476	act = 'READ'
477	pos = 'REWIND'
478	! if( form_flag.EQ.MPP_NETCDF )call mpp_error( FATAL, 'MPP_OPEN: only writes are currently supported with netCDF.' )
479	else if( action_flag.EQ.MPP_WRONLY .OR. action_flag.EQ.MPP_OVERWR )then
480	act = 'WRITE'
481	pos = 'REWIND'
482	else if( action_flag.EQ.MPP_APPEND )then
483	act = 'WRITE'
484	pos = 'APPEND'
485	else
486	call mpp_error( FATAL, 'MPP_OPEN: action must be one of MPP_WRONLY, MPP_APPEND or MPP_RDONLY.' )
487	end if
488
489	!access: sequential or direct: ignored by netCDF
490	if( form_flag.NE.MPP_NETCDF )then
491	if( access_flag.EQ.MPP_SEQUENTIAL )then
492	acc = 'SEQUENTIAL'
493	else if( access_flag.EQ.MPP_DIRECT )then
494	acc = 'DIRECT'
495	if( form_flag.EQ.MPP_ASCII )call mpp_error( FATAL, 'MPP_OPEN: formatted direct access I/O is prohibited.' )
496	if( .NOT.PRESENT(recl) ) &
497	call mpp_error( FATAL, 'MPP_OPEN: recl (record length in bytes) must be specified with access=MPP_DIRECT.' )
498	mpp_file(unit)%record = 1
499	records_per_pe = 1 !each PE writes 1 record per mpp_write
500	else
501	call mpp_error( FATAL, 'MPP_OPEN: access must be one of MPP_SEQUENTIAL or MPP_DIRECT.' )
502	end if
503	end if
504
505	!threading: SINGLE or MULTI
506	if( threading_flag.EQ.MPP_MULTI )then
507	!fileset: MULTI or SINGLE (only for multi-threaded I/O
508	if( fileset_flag.EQ.MPP_SINGLE )then
509	if( form_flag.EQ.MPP_NETCDF .AND. act.EQ.'WRITE' ) &
510	call mpp_error( FATAL, 'MPP_OPEN: netCDF currently does not support single-file multi-threaded output.' )
511
512	#ifdef _CRAYT3E
513	call ASSIGN( 'assign -I -F global.privpos f:'//trim(mpp_file(unit)%name), error )
514	#endif
515	else if( fileset_flag.NE.MPP_MULTI )then
516	call mpp_error( FATAL, 'MPP_OPEN: fileset must be one of MPP_MULTI or MPP_SINGLE.' )
517	end if
518	else if( threading_flag.NE.MPP_SINGLE )then
519	call mpp_error( FATAL, 'MPP_OPEN: threading must be one of MPP_SINGLE or MPP_MULTI.' )
520	end if
521
522	!apply I/O specs before opening the file
523	!note that -P refers to the scope of a fortran unit, which is always thread-private even if file is shared
524	#ifdef CRAYPVP
525	call ASSIGN( 'assign -I -P thread f:'//trim(mpp_file(unit)%name), error )
526	#endif
527	#ifdef _CRAYT3E
528	call ASSIGN( 'assign -I -P private f:'//trim(mpp_file(unit)%name), error )
529	#endif
530	if( PRESENT(iospec) )then
531	!iospec provides hints to the system on how to organize I/O
532	!on Cray systems this is done through 'assign', see assign(1) and assign(3F)
533	!on other systems this will be expanded as needed
534	!no error checks here on whether the supplied iospec is valid
535	#ifdef SGICRAY
536	call ASSIGN( 'assign -I '//trim(iospec)//' f:'//trim(mpp_file(unit)%name), error )
537	if( form_flag.EQ.MPP_NETCDF )then
538	!for netCDF on SGI/Cray systems we pass it to the environment variable NETCDF_XFFIOSPEC
539	!ideally we should parse iospec, pass the argument of -F to NETCDF_FFIOSPEC, and the rest to NETCDF_XFFIOSPEC
540	!maybe I'll get around to it someday
541	!PXFSETENV is a POSIX-standard routine for setting environment variables from fortran
542	call PXFSETENV( 'NETCDF_XFFIOSPEC', 0, trim(iospec), 0, 1, error )
543	end if
544	#endif
545	end if
546
547	!open the file as specified above for various formats
548	if( form_flag.EQ.MPP_NETCDF )then
549	#ifdef use_netCDF
550	if( action_flag.EQ.MPP_WRONLY )then
551	error = NF_CREATE( trim(mpp_file(unit)%name), NF_NOCLOBBER, mpp_file(unit)%ncid ); call netcdf_err(error)
552	if( verbose )print '(a,i3,i16)', 'MPP_OPEN: new netCDF file: pe, ncid=', pe, mpp_file(unit)%ncid
553	else if( action_flag.EQ.MPP_OVERWR )then
554	error = NF_CREATE( trim(mpp_file(unit)%name), NF_CLOBBER, mpp_file(unit)%ncid ); call netcdf_err(error)
555	action_flag = MPP_WRONLY !after setting clobber, there is no further distinction btwn MPP_WRONLY and MPP_OVERWR
556	if( verbose )print '(a,i3,i16)', 'MPP_OPEN: overwrite netCDF file: pe, ncid=', pe, mpp_file(unit)%ncid
557	else if( action_flag.EQ.MPP_APPEND )then
558	error = NF_OPEN( trim(mpp_file(unit)%name), NF_WRITE, mpp_file(unit)%ncid ); call netcdf_err(error)
559	!get the current time level of the file: writes to this file will be at next time level
560	error = NF_INQ_UNLIMDIM( mpp_file(unit)%ncid, unlim%did )
561	if( error.EQ.NF_NOERR )then
562	error = NF_INQ_DIM( mpp_file(unit)%ncid, unlim%did, unlim%name, mpp_file(unit)%time_level )
563	call netcdf_err(error)
564	error = NF_INQ_VARID( mpp_file(unit)%ncid, unlim%name, mpp_file(unit)%id ); call netcdf_err(error)
565	end if
566	if( verbose )print '(a,i3,i16,i4)', 'MPP_OPEN: append to existing netCDF file: pe, ncid, time_axis_id=',&
567	pe, mpp_file(unit)%ncid, mpp_file(unit)%id
568	else if( action_flag.EQ.MPP_RDONLY )then
569	error = NF_OPEN( trim(mpp_file(unit)%name), NF_NOWRITE, mpp_file(unit)%ncid ); call netcdf_err(error)
570	if( verbose )print '(a,i3,i16,i4)', 'MPP_OPEN: opening existing netCDF file: pe, ncid, time_axis_id=',&
571	pe, mpp_file(unit)%ncid, mpp_file(unit)%id
572	mpp_file(unit)%format=form_flag ! need this for mpp_read
573	call mpp_read_meta(unit)
574	end if
575	mpp_file(unit)%opened = .TRUE.
576	#endif
577	else
578	!format: ascii, native, or IEEE 32 bit
579	if( form_flag.EQ.MPP_ASCII )then
580	for = 'FORMATTED'
581	else if( form_flag.EQ.MPP_IEEE32 )then
582	for = 'UNFORMATTED'
583	!assign -N is currently unsupported on SGI
584	#ifdef _CRAY
585	call ASSIGN( 'assign -I -N ieee_32 f:'//trim(mpp_file(unit)%name), error )
586	#endif
587	else if( form_flag.EQ.MPP_NATIVE )then
588	for = 'UNFORMATTED'
589	else
590	call mpp_error( FATAL, 'MPP_OPEN: form must be one of MPP_ASCII, MPP_NATIVE, MPP_IEEE32 or MPP_NETCDF.' )
591	end if
592	inquire( file=trim(mpp_file(unit)%name), EXIST=exists )
593	if( exists .AND. action_flag.EQ.MPP_WRONLY ) &
594	call mpp_error( WARNING, 'MPP_OPEN: File '//trim(mpp_file(unit)%name)//' opened WRONLY already exists!' )
595	if( action_flag.EQ.MPP_OVERWR )action_flag = MPP_WRONLY
596	!perform the OPEN here
597	if( PRESENT(recl) )then
598	if( verbose )print '(2(x,a,i3),5(x,a),a,i8)', 'MPP_OPEN: PE=', pe, &
599	'unit=', unit, trim(mpp_file(unit)%name), 'attributes=', trim(acc), trim(for), trim(act), ' RECL=', recl
600	open( unit, file=trim(mpp_file(unit)%name), access=acc, form=for, action=act, recl=recl )
601	else
602	if( verbose )print '(2(x,a,i3),6(x,a))', 'MPP_OPEN: PE=', pe, &
603	'unit=', unit, trim(mpp_file(unit)%name), 'attributes=', trim(acc), trim(for), trim(pos), trim(act)
604	open( unit, file=trim(mpp_file(unit)%name), access=acc, form=for, action=act, position=pos )
605	end if
606	!check if OPEN worked
607	inquire( unit,OPENED=mpp_file(unit)%opened )
608	if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_OPEN: error in OPEN() statement.' )
609	end if
610	mpp_file(unit)%action = action_flag
611	mpp_file(unit)%format = form_flag
612	mpp_file(unit)%access = access_flag
613	mpp_file(unit)%threading = threading_flag
614	mpp_file(unit)%fileset = fileset_flag
615	if( PRESENT(nohdrs) )mpp_file(unit)%nohdrs = nohdrs
616
617	if( action_flag.EQ.MPP_WRONLY )then
618	if( form_flag.NE.MPP_NETCDF .AND. access_flag.EQ.MPP_DIRECT )call mpp_write_meta( unit, 'record_length', ival=recl )
619	!actual file name
620	call mpp_write_meta( unit, 'filename', cval=mpp_file(unit)%name )
621	!MPP_IO package version
622	call mpp_write_meta( unit, 'MPP_IO_VERSION', cval=trim(version) )
623	!filecount for multifileset
624	if( threading_flag.EQ.MPP_MULTI .AND. fileset_flag.EQ.MPP_MULTI ) &
625	call mpp_write_meta( unit, 'NumFilesInSet', ival=npes )
626	end if
627
628	return
629	end subroutine mpp_open
630
631	subroutine mpp_close( unit, action )
632	integer, intent(in) :: unit
633	integer, intent(in), optional :: action
634	character(len=8) :: status
635	logical :: collect
636
637	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_CLOSE: must first call mpp_io_init.' )
638	if( unit.EQ.NULLUNIT )return !nothing was actually opened on this unit
639
640	!action on close
641	status = 'KEEP'
642	!collect is supposed to launch the post-processing collector tool for multi-fileset
643	collect = .FALSE.
644	if( PRESENT(action) )then
645	if( action.EQ.MPP_DELETE )then
646	status = 'DELETE'
647	else if( action.EQ.MPP_COLLECT )then
648	collect = .FALSE. !should be TRUE but this is not yet ready
649	call mpp_error( WARNING, 'MPP_CLOSE: the COLLECT operation is not yet implemented.' )
650	else
651	call mpp_error( FATAL, 'MPP_CLOSE: action must be one of MPP_DELETE or MPP_COLLECT.' )
652	end if
653	end if
654	if( mpp_file(unit)%fileset.NE.MPP_MULTI )collect = .FALSE.
655	if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
656	#ifdef use_netCDF
657	error = NF_CLOSE(mpp_file(unit)%ncid); call netcdf_err(error)
658	#endif
659	else
660	close(unit,status=status)
661	end if
662	#ifdef SGICRAY
663	!this line deleted: since the FILENV is a shared file, this might cause a problem in
664	! multi-threaded I/O if one PE does assign -R before another one has opened it.
665	! call ASSIGN( 'assign -R f:'//trim(mpp_file(unit)%name), error )
666	#endif
667	mpp_file(unit)%name = ' '
668	mpp_file(unit)%action = -1
669	mpp_file(unit)%format = -1
670	mpp_file(unit)%access = -1
671	mpp_file(unit)%threading = -1
672	mpp_file(unit)%fileset = -1
673	mpp_file(unit)%record = -1
674	mpp_file(unit)%ncid = -1
675	mpp_file(unit)%opened = .FALSE.
676	mpp_file(unit)%initialized = .FALSE.
677	mpp_file(unit)%id = -1
678	mpp_file(unit)%time_level = 0
679	mpp_file(unit)%time = NULLTIME
680	return
681	end subroutine mpp_close
682
683	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
684	! !
685	! MPP_WRITE_META !
686	! !
687	!This series of routines is used to describe the contents of the file !
688	!being written on <unit>. Each file can contain any number of fields, !
689	!which can be functions of 0-3 spatial axes and 0-1 time axes. Axis !
690	!descriptors are stored in the <axistype> structure and field !
691	!descriptors in the <fieldtype> structure. !
692	! !
693	! type, public :: axistype !
694	! sequence !
695	! character(len=128) :: name !
696	! character(len=128) :: units !
697	! character(len=256) :: longname !
698	! integer :: sense !+/-1, depth or height? !
699	! type(domain1D) :: domain !
700	! real, pointer :: data(:) !axis values (not used if time axis) !
701	! integer :: id !
702	! end type axistype !
703	! !
704	! type, public :: fieldtype !
705	! sequence !
706	! character(len=128) :: name !
707	! character(len=128) :: units !
708	! character(len=256) :: longname !
709	! real :: min, max, missing, fill, scale, add !
710	! type(axistype), pointer :: axis(:) !
711	! integer :: id !
712	! end type fieldtype !
713	! !
714	!The metadata contained in the type is always written for each axis and !
715	!field. Any other metadata one wishes to attach to an axis or field !
716	!can subsequently be passed to mpp_write_meta using the ID, as shown below. !
717	! !
718	!mpp_write_meta can take several forms: !
719	! !
720	! mpp_write_meta( unit, name, rval=rval, pack=pack ) !
721	! mpp_write_meta( unit, name, ival=ival ) !
722	! mpp_write_meta( unit, name, cval=cval ) !
723	! integer, intent(in) :: unit !
724	! character(len=*), intent(in) :: name !
725	! real, intent(in), optional :: rval(:) !
726	! integer, intent(in), optional :: ival(:) !
727	! character(len=*), intent(in), optional :: cval !
728	! !
729	! This form defines global metadata associated with the file as a !
730	! whole. The attribute is named <name> and can take on a real, integer !
731	! or character value. <rval> and <ival> can be scalar or 1D arrays. !
732	! !
733	! mpp_write_meta( unit, id, name, rval=rval, pack=pack ) !
734	! mpp_write_meta( unit, id, name, ival=ival ) !
735	! mpp_write_meta( unit, id, name, cval=cval ) !
736	! integer, intent(in) :: unit, id !
737	! character(len=*), intent(in) :: name !
738	! real, intent(in), optional :: rval(:) !
739	! integer, intent(in), optional :: ival(:) !
740	! character(len=*), intent(in), optional :: cval !
741	! !
742	! This form defines metadata associated with a previously defined !
743	! axis or field, identified to mpp_write_meta by its unique ID <id>. !
744	! The attribute is named <name> and can take on a real, integer !
745	! or character value. <rval> and <ival> can be scalar or 1D arrays. !
746	! This need not be called for attributes already contained in !
747	! the type. !
748	! !
749	! PACK can take values 1,2,4,8. This only has meaning when writing !
750	! floating point numbers. The value of PACK defines the number of words !
751	! written into 8 bytes. For pack=4 and pack=8, an integer value is !
752	! written: rval is assumed to have been scaled to the appropriate dynamic !
753	! range. !
754	! PACK currently only works for netCDF files, and is ignored otherwise. !
755	! !
756	! subroutine mpp_write_meta_axis( unit, axis, name, units, longname, & !
757	! cartesian, sense, domain, data ) !
758	! integer, intent(in) :: unit !
759	! type(axistype), intent(inout) :: axis !
760	! character(len=*), intent(in) :: name, units, longname !
761	! character(len=*), intent(in), optional :: cartesian !
762	! integer, intent(in), optional :: sense !
763	! type(domain1D), intent(in), optional :: domain !
764	! real, intent(in), optional :: data(:) !
765	! !
766	! This form defines a time or space axis. Metadata corresponding to the !
767	! type above are written to the file on <unit>. A unique ID for subsequent!
768	! references to this axis is returned in axis%id. If the <domain> !
769	! element is present, this is recognized as a distributed data axis !
770	! and domain decomposition information is also written if required (the !
771	! domain decomposition info is required for multi-fileset multi-threaded !
772	! I/O). If the <data> element is allocated, it is considered to be a space!
773	! axis, otherwise it is a time axis with an unlimited dimension. Only one !
774	! time axis is allowed per file. !
775	! !
776	! subroutine mpp_write_meta_field( unit, field, axes, name, units, longname!
777	! min, max, missing, fill, scale, add, pack ) !
778	! integer, intent(in) :: unit !
779	! type(fieldtype), intent(out) :: field !
780	! type(axistype), intent(in) :: axes(:) !
781	! character(len=*), intent(in) :: name, units, longname !
782	! real, intent(in), optional :: min, max, missing, fill, scale, add !
783	! integer, intent(in), optional :: pack !
784	! !
785	! This form defines a field. Metadata corresponding to the type !
786	! above are written to the file on <unit>. A unique ID for subsequent !
787	! references to this field is returned in field%id. At least one axis !
788	! must be associated, 0D variables are not considered. mpp_write_meta !
789	! must previously have been called on all axes associated with this !
790	! field. !
791	! !
792	! The mpp_write_meta package also includes subroutines write_attribute and !
793	! write_attribute_netcdf, that are private to this module. !
794	! !
795	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
796	subroutine mpp_write_meta_global( unit, name, rval, ival, cval, pack )
797	!writes a global metadata attribute to unit <unit>
798	!attribute <name> can be an real, integer or character
799	!one and only one of rval, ival, and cval should be present
800	!the first found will be used
801	!for a non-netCDF file, it is encoded into a string "GLOBAL <name> <val>"
802	integer, intent(in) :: unit
803	character(len=*), intent(in) :: name
804	real, intent(in), optional :: rval(:)
805	integer, intent(in), optional :: ival(:)
806	character(len=*), intent(in), optional :: cval
807	integer, intent(in), optional :: pack
808
809	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_WRITE_META: must first call mpp_io_init.' )
810	if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_WRITE_META: invalid unit number.' )
811	if( mpp_file(unit)%threading.EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
812	if( mpp_file(unit)%fileset.EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
813	if( mpp_file(unit)%action.NE.MPP_WRONLY )return !no writing metadata on APPEND
814	if( mpp_file(unit)%initialized ) &
815	call mpp_error( FATAL, 'MPP_WRITE_META: cannot write metadata to file after an mpp_write.' )
816
817	if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
818	#ifdef use_netCDF
819	call write_attribute_netcdf( unit, NF_GLOBAL, name, rval, ival, cval, pack )
820	#endif
821	else
822	call write_attribute( unit, 'GLOBAL '//trim(name), rval, ival, cval, pack )
823	end if
824
825	return
826	end subroutine mpp_write_meta_global
827
828	!versions of above to support <rval> and <ival> as scalars (because of f90 strict rank matching)
829	subroutine mpp_write_meta_global_scalar_r( unit, name, rval, pack )
830	integer, intent(in) :: unit
831	character(len=*), intent(in) :: name
832	real, intent(in) :: rval
833	integer, intent(in), optional :: pack
834
835	call mpp_write_meta_global( unit, name, rval=(/rval/), pack=pack )
836	return
837	end subroutine mpp_write_meta_global_scalar_r
838
839	subroutine mpp_write_meta_global_scalar_i( unit, name, ival )
840	integer, intent(in) :: unit
841	character(len=*), intent(in) :: name
842	integer, intent(in) :: ival
843
844	call mpp_write_meta_global( unit, name, ival=(/ival/) )
845	return
846	end subroutine mpp_write_meta_global_scalar_i
847
848	subroutine mpp_write_meta_var( unit, id, name, rval, ival, cval, pack )
849	!writes a metadata attribute for variable <id> to unit <unit>
850	!attribute <name> can be an real, integer or character
851	!one and only one of rval, ival, and cval should be present
852	!the first found will be used
853	!for a non-netCDF file, it is encoded into a string "<id> <name> <val>"
854	integer, intent(in) :: unit, id
855	character(len=*), intent(in) :: name
856	real, intent(in), optional :: rval(:)
857	integer, intent(in), optional :: ival(:)
858	character(len=*), intent(in), optional :: cval
859	integer, intent(in), optional :: pack
860
861	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_WRITE_META: must first call mpp_io_init.' )
862	if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_WRITE_META: invalid unit number.' )
863	if( mpp_file(unit)%threading.EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
864	if( mpp_file(unit)%fileset.EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
865	if( mpp_file(unit)%action.NE.MPP_WRONLY )return !no writing metadata on APPEND
866	if( mpp_file(unit)%initialized ) &
867	call mpp_error( FATAL, 'MPP_WRITE_META: cannot write metadata to file after an mpp_write.' )
868
869	if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
870	call write_attribute_netcdf( unit, id, name, rval, ival, cval, pack )
871	else
872	write( text, '(a,i4,a)' )'VARIABLE ', id, ' '//name
873	call write_attribute( unit, trim(text), rval, ival, cval, pack )
874	end if
875
876	return
877	end subroutine mpp_write_meta_var
878
879	!versions of above to support <rval> and <ival> as scalar (because of f90 strict rank matching)
880	subroutine mpp_write_meta_scalar_r( unit, id, name, rval, pack )
881	integer, intent(in) :: unit, id
882	character(len=*), intent(in) :: name
883	real, intent(in) :: rval
884	integer, intent(in), optional :: pack
885
886	call mpp_write_meta( unit, id, name, rval=(/rval/), pack=pack )
887	return
888	end subroutine mpp_write_meta_scalar_r
889
890	subroutine mpp_write_meta_scalar_i( unit, id, name, ival )
891	integer, intent(in) :: unit, id
892	character(len=*), intent(in) :: name
893	integer, intent(in) :: ival
894
895	call mpp_write_meta( unit, id, name, ival=(/ival/) )
896	return
897	end subroutine mpp_write_meta_scalar_i
898
899	subroutine mpp_write_meta_axis( unit, axis, name, units, longname, cartesian, sense, domain, data )
900	!load the values in an axistype (still need to call mpp_write)
901	!write metadata attributes for axis
902	!it is declared intent(inout) so you can nullify pointers in the incoming object if needed
903	!the f90 standard doesn't guarantee that intent(out) on a type guarantees that its pointer components will be unassociated
904	integer, intent(in) :: unit
905	type(axistype), intent(inout) :: axis
906	character(len=*), intent(in) :: name, units, longname
907	character(len=*), intent(in), optional :: cartesian
908	integer, intent(in), optional :: sense
909	type(domain1D), intent(in), optional :: domain
910	real, intent(in), optional :: data(:)
911	integer :: is, ie, isg, ieg
912
913	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_WRITE_META: must first call mpp_io_init.' )
914	if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_WRITE_META: invalid unit number.' )
915	if( mpp_file(unit)%threading.EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
916	if( mpp_file(unit)%fileset.EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
917	if( mpp_file(unit)%action.NE.MPP_WRONLY )return !no writing metadata on APPEND
918	if( mpp_file(unit)%initialized ) &
919	call mpp_error( FATAL, 'MPP_WRITE_META: cannot write metadata to file after an mpp_write.' )
920
921	!pre-existing pointers need to be nullified
922	if( ASSOCIATED(axis%data) )NULLIFY(axis%data)
923	!load axistype
924	axis%name = name
925	axis%units = units
926	axis%longname = longname
927	if( PRESENT(cartesian) )axis%cartesian = cartesian
928	if( PRESENT(sense) )axis%sense = sense
929	if( PRESENT(domain) )then
930	axis%domain = domain
931	call mpp_get_global_domain( domain, isg, ieg )
932	call mpp_get_compute_domain( domain, is, ie )
933	else
934	axis%domain = NULL_DOMAIN1D
935	if( PRESENT(data) )then
936	isg=1; ieg=size(data); is=isg; ie=ieg
937	endif
938	end if
939	if( PRESENT(data) )then
940	if( PRESENT(domain) )then
941	if( size(data).NE.ieg-isg+1 ) &
942	call mpp_error( FATAL, 'MPP_WRITE_META_AXIS: size(data).NE.domain%global%size.' )
943	allocate( axis%data(isg:ieg) )
944	else
945	allocate( axis%data(size(data)) )
946	end if
947	axis%data = data
948	end if
949
950	!write metadata
951	if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
952	#ifdef use_netCDF
953	!write axis def
954	!space axes are always floats, time axis is always double
955	if( ASSOCIATED(axis%data) )then !space axis
956	if( mpp_file(unit)%fileset.EQ.MPP_MULTI .AND. axis%domain.NE.NULL_DOMAIN1D )then
957	error = NF_DEF_DIM( mpp_file(unit)%ncid, axis%name, ie-is+1, axis%did )
958	else
959	error = NF_DEF_DIM( mpp_file(unit)%ncid, axis%name, size(axis%data), axis%did )
960	end if
961	call netcdf_err(error)
962	error = NF_DEF_VAR( mpp_file(unit)%ncid, axis%name, NF_FLOAT, 1, axis%did, axis%id ); call netcdf_err(error)
963	else !time axis
964	if( mpp_file(unit)%id.NE.-1 ) &
965	call mpp_error( FATAL, 'MPP_WRITE_META_AXIS: There is already a time axis for this file.' )
966	error = NF_DEF_DIM( mpp_file(unit)%ncid, axis%name, NF_UNLIMITED, axis%did ); call netcdf_err(error)
967	error = NF_DEF_VAR( mpp_file(unit)%ncid, axis%name, NF_DOUBLE, 1, axis%did, axis%id ); call netcdf_err(error)
968	mpp_file(unit)%id = axis%id !file ID is the same as time axis varID
969	end if
970	#endif
971	else
972	varnum = varnum + 1
973	axis%id = varnum
974	axis%did = varnum
975	!write axis def
976	write( text, '(a,i4,a)' )'AXIS ', axis%id, ' name'
977	call write_attribute( unit, trim(text), cval=axis%name )
978	write( text, '(a,i4,a)' )'AXIS ', axis%id, ' size'
979	if( ASSOCIATED(axis%data) )then !space axis
980	if( mpp_file(unit)%fileset.EQ.MPP_MULTI .AND. axis%domain.NE.NULL_DOMAIN1D )then
981	call write_attribute( unit, trim(text), ival=(/ie-is+1/) )
982	else
983	call write_attribute( unit, trim(text), ival=(/size(axis%data)/) )
984	end if
985	else !time axis
986	if( mpp_file(unit)%id.NE.-1 ) &
987	call mpp_error( FATAL, 'MPP_WRITE_META_AXIS: There is already a time axis for this file.' )
988	call write_attribute( unit, trim(text), ival=(/0/) ) !a size of 0 indicates time axis
989	mpp_file(unit)%id = axis%id
990	end if
991	end if
992	!write axis attributes
993	call mpp_write_meta( unit, axis%id, 'long_name', cval=axis%longname )
994	call mpp_write_meta( unit, axis%id, 'units', cval=axis%units )
995	if( PRESENT(cartesian) )call mpp_write_meta( unit, axis%id, 'cartesian_axis', cval=axis%cartesian )
996	if( PRESENT(sense) )then
997	if( sense.EQ.-1 )then
998	call mpp_write_meta( unit, axis%id, 'positive', cval='down' )
999	else if( sense.EQ.1 )then
1000	call mpp_write_meta( unit, axis%id, 'positive', cval='up' )
1001	end if
1002	!silently ignore values of sense other than +/-1.
1003	end if
1004	if( mpp_file(unit)%threading.EQ.MPP_MULTI .AND. mpp_file(unit)%fileset.EQ.MPP_MULTI .AND. axis%domain.NE.NULL_DOMAIN1D )then
1005	call mpp_write_meta( unit, axis%id, 'domain_decomposition', ival=(/isg,ieg,is,ie/) )
1006	end if
1007	if( verbose )print '(a,2i3,x,a,2i3)', 'MPP_WRITE_META: Wrote axis metadata, pe, unit, axis%name, axis%id, axis%did=', &
1008	pe, unit, trim(axis%name), axis%id, axis%did
1009
1010	return
1011	end subroutine mpp_write_meta_axis
1012
1013	subroutine mpp_write_meta_field( unit, field, axes, name, units, longname, min, max, missing, fill, scale, add, pack )
1014	!define field: must have already called mpp_write_meta(axis) for each axis
1015	integer, intent(in) :: unit
1016	type(fieldtype), intent(out) :: field
1017	type(axistype), intent(in) :: axes(:)
1018	character(len=*), intent(in) :: name, units, longname
1019	real, intent(in), optional :: min, max, missing, fill, scale, add
1020	integer, intent(in), optional :: pack
1021	!this array is required because of f77 binding on netCDF interface
1022	integer, allocatable :: axis_id(:)
1023	real :: a, b
1024	integer :: i
1025
1026	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_WRITE_META: must first call mpp_io_init.' )
1027	if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_WRITE_META: invalid unit number.' )
1028	if( mpp_file(unit)%threading.EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
1029	if( mpp_file(unit)%fileset.EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
1030	if( mpp_file(unit)%action.NE.MPP_WRONLY )return !no writing metadata on APPEND
1031	if( mpp_file(unit)%initialized ) &
1032	call mpp_error( FATAL, 'MPP_WRITE_META: cannot write metadata to file after an mpp_write.' )
1033
1034	!pre-existing pointers need to be nullified
1035	if( ASSOCIATED(field%axes) )NULLIFY(field%axes)
1036	!fill in field metadata
1037	field%name = name
1038	field%units = units
1039	field%longname = longname
1040	allocate( field%axes(size(axes)) )
1041	field%axes = axes
1042	field%time_axis_index = -1 !this value will never match any axis index
1043	!size is buffer area for the corresponding axis info: it is required to buffer this info in the fieldtype
1044	!because axis might be reused in different files
1045	allocate( field%size(size(axes)) )
1046	do i = 1,size(axes)
1047	if( ASSOCIATED(axes(i)%data) )then !space axis
1048	field%size(i) = size(axes(i)%data)
1049	else !time
1050	field%size(i) = 1
1051	field%time_axis_index = i
1052	end if
1053	end do
1054	!attributes
1055	if( PRESENT(min) )field%min = min
1056	if( PRESENT(max) )field%max = max
1057	if( PRESENT(missing) )field%missing = missing
1058	if( PRESENT(fill) )field%fill = fill
1059	if( PRESENT(scale) )field%scale = scale
1060	if( PRESENT(add) )field%add = add
1061
1062	!pack is currently used only for netCDF
1063	field%pack = 2 !default write 32-bit floats
1064	if( PRESENT(pack) )field%pack = pack
1065	if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
1066	#ifdef use_netCDF
1067	allocate( axis_id(size(field%axes)) )
1068	do i = 1,size(field%axes)
1069	axis_id(i) = field%axes(i)%did
1070	end do
1071	!write field def
1072	select case (field%pack)
1073	case(1)
1074	error = NF_DEF_VAR( mpp_file(unit)%ncid, field%name, NF_DOUBLE, size(field%axes), axis_id, field%id )
1075	case(2)
1076	error = NF_DEF_VAR( mpp_file(unit)%ncid, field%name, NF_FLOAT, size(field%axes), axis_id, field%id )
1077	case(4)
1078	if( .NOT.PRESENT(scale) .OR. .NOT.PRESENT(add) ) &
1079	call mpp_error( FATAL, 'MPP_WRITE_META_FIELD: scale and add must be supplied when pack=4.' )
1080	error = NF_DEF_VAR( mpp_file(unit)%ncid, field%name, NF_SHORT, size(field%axes), axis_id, field%id )
1081	case(8)
1082	if( .NOT.PRESENT(scale) .OR. .NOT.PRESENT(add) ) &
1083	call mpp_error( FATAL, 'MPP_WRITE_META_FIELD: scale and add must be supplied when pack=8.' )
1084	error = NF_DEF_VAR( mpp_file(unit)%ncid, field%name, NF_BYTE, size(field%axes), axis_id, field%id )
1085	case default
1086	call mpp_error( FATAL, 'MPP_WRITE_META_FIELD: only legal packing values are 1,2,4,8.' )
1087	end select
1088	call netcdf_err(error)
1089	#endif
1090	else
1091	varnum = varnum + 1
1092	field%id = varnum
1093	if( PRESENT(pack) )call mpp_error( WARNING, 'MPP_WRITE_META: Packing is currently available only on netCDF files.' )
1094	!write field def
1095	write( text, '(a,i4,a)' )'FIELD ', field%id, ' name'
1096	call write_attribute( unit, trim(text), cval=field%name )
1097	write( text, '(a,i4,a)' )'FIELD ', field%id, ' axes'
1098	call write_attribute( unit, trim(text), ival=field%axes(:)%did )
1099	end if
1100	!write field attributes: these names follow netCDF conventions
1101	call mpp_write_meta( unit, field%id, 'long_name', cval=field%longname )
1102	call mpp_write_meta( unit, field%id, 'units', cval=field%units )
1103	!all real attributes must be written as packed
1104	if( PRESENT(min) .AND. PRESENT(max) )then
1105	if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
1106	call mpp_write_meta( unit, field%id, 'valid_range', rval=(/min,max/), pack=pack )
1107	else
1108	a = nint((min-add)/scale)
1109	b = nint((max-add)/scale)
1110	call mpp_write_meta( unit, field%id, 'valid_range', rval=(/a, b /), pack=pack )
1111	end if
1112	else if( PRESENT(min) )then
1113	if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
1114	call mpp_write_meta( unit, field%id, 'valid_min', rval=field%min, pack=pack )
1115	else
1116	a = nint((min-add)/scale)
1117	call mpp_write_meta( unit, field%id, 'valid_min', rval=a, pack=pack )
1118	end if
1119	else if( PRESENT(max) )then
1120	if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
1121	call mpp_write_meta( unit, field%id, 'valid_max', rval=field%max, pack=pack )
1122	else
1123	a = nint((max-add)/scale)
1124	call mpp_write_meta( unit, field%id, 'valid_max', rval=a, pack=pack )
1125	end if
1126	end if
1127	if( PRESENT(missing) )then
1128	if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
1129	call mpp_write_meta( unit, field%id, 'missing_value', rval=field%missing, pack=pack )
1130	else
1131	a = nint((missing-add)/scale)
1132	call mpp_write_meta( unit, field%id, 'missing_value', rval=a, pack=pack )
1133	end if
1134	end if
1135	if( PRESENT(fill) )then
1136	if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
1137	call mpp_write_meta( unit, field%id, '_FillValue', rval=field%missing, pack=pack )
1138	else
1139	a = nint((fill-add)/scale)
1140	call mpp_write_meta( unit, field%id, '_FillValue', rval=a, pack=pack )
1141	end if
1142	end if
1143	if( field%pack.NE.1 .AND. field%pack.NE.2 )then
1144	call mpp_write_meta( unit, field%id, 'packing', ival=field%pack )
1145	if( PRESENT(scale) )call mpp_write_meta( unit, field%id, 'scale_factor', rval=field%scale )
1146	if( PRESENT(add) )call mpp_write_meta( unit, field%id, 'add_offset', rval=field%add )
1147	end if
1148	if( verbose )print '(a,2i3,x,a,i3)', 'MPP_WRITE_META: Wrote field metadata: pe, unit, field%name, field%id=', &
1149	pe, unit, trim(field%name), field%id
1150
1151	return
1152	end subroutine mpp_write_meta_field
1153
1154	subroutine write_attribute( unit, name, rval, ival, cval, pack )
1155	!called to write metadata for non-netCDF I/O
1156	integer, intent(in) :: unit
1157	character(len=*), intent(in) :: name
1158	real, intent(in), optional :: rval(:)
1159	integer, intent(in), optional :: ival(:)
1160	character(len=*), intent(in), optional :: cval
1161	!pack is currently ignored in this routine: only used by netCDF I/O
1162	integer, intent(in), optional :: pack
1163
1164	if( mpp_file(unit)%nohdrs )return
1165	!encode text string
1166	if( PRESENT(rval) )then
1167	write( text,* )trim(name)//'=', rval
1168	else if( PRESENT(ival) )then
1169	write( text,* )trim(name)//'=', ival
1170	else if( PRESENT(cval) )then
1171	text = ' '//trim(name)//'='//trim(cval)
1172	else
1173	call mpp_error( FATAL, 'WRITE_ATTRIBUTE: one of rval, ival, cval must be present.' )
1174	end if
1175	if( mpp_file(unit)%format.EQ.MPP_ASCII )then
1176	!implies sequential access
1177	write( unit,fmt='(a)' )trim(text)//char(10)
1178	else !MPP_IEEE32 or MPP_NATIVE
1179	if( mpp_file(unit)%access.EQ.MPP_SEQUENTIAL )then
1180	write(unit)trim(text)//char(10)
1181	else !MPP_DIRECT
1182	write( unit,rec=mpp_file(unit)%record )trim(text)//char(10)
1183	if( verbose )print '(a,i3,a,i3)', 'WRITE_ATTRIBUTE: PE=', pe, ' wrote record ', mpp_file(unit)%record
1184	mpp_file(unit)%record = mpp_file(unit)%record + 1
1185	end if
1186	end if
1187	return
1188	end subroutine write_attribute
1189
1190	subroutine write_attribute_netcdf( unit, id, name, rval, ival, cval, pack )
1191	!called to write metadata for netCDF I/O
1192	integer, intent(in) :: unit
1193	integer, intent(in) :: id
1194	character(len=*), intent(in) :: name
1195	real, intent(in), optional :: rval(:)
1196	integer, intent(in), optional :: ival(:)
1197	character(len=*), intent(in), optional :: cval
1198	integer, intent(in), optional :: pack
1199	integer :: lenc
1200	integer, allocatable :: rval_i(:)
1201	#ifdef use_netCDF
1202	if( PRESENT(rval) )then
1203	!pack is only meaningful for FP numbers
1204	if( PRESENT(pack) )then
1205	if( pack.EQ.1 )then
1206	if( KIND(rval).EQ.DOUBLE_KIND )then
1207	error = NF_PUT_ATT_DOUBLE( mpp_file(unit)%ncid, id, name, NF_DOUBLE, size(rval), rval )
1208	else if( KIND(rval).EQ.FLOAT_KIND )then
1209	call mpp_error( WARNING, &
1210	'WRITE_ATTRIBUTE_NETCDF: attempting to write internal 32-bit real as external 64-bit.' )
1211	error = NF_PUT_ATT_REAL ( mpp_file(unit)%ncid, id, name, NF_DOUBLE, size(rval), rval )
1212	end if
1213	call netcdf_err(error)
1214	else if( pack.EQ.2 )then
1215	if( KIND(rval).EQ.DOUBLE_KIND )then
1216	error = NF_PUT_ATT_DOUBLE( mpp_file(unit)%ncid, id, name, NF_FLOAT, size(rval), rval )
1217	else if( KIND(rval).EQ.FLOAT_KIND )then
1218	error = NF_PUT_ATT_REAL ( mpp_file(unit)%ncid, id, name, NF_FLOAT, size(rval), rval )
1219	end if
1220	call netcdf_err(error)
1221	else if( pack.EQ.4 )then
1222	allocate( rval_i(size(rval)) )
1223	rval_i = rval
1224	if( KIND(rval).EQ.DOUBLE_KIND )then
1225	error = NF_PUT_ATT_DOUBLE( mpp_file(unit)%ncid, id, name, NF_SHORT, size(rval_i), rval )
1226	else if( KIND(rval).EQ.FLOAT_KIND )then
1227	error = NF_PUT_ATT_REAL ( mpp_file(unit)%ncid, id, name, NF_SHORT, size(rval_i), rval )
1228	end if
1229	call netcdf_err(error)
1230	deallocate(rval_i)
1231	else if( pack.EQ.8 )then
1232	allocate( rval_i(size(rval)) )
1233	rval_i = rval
1234	if( KIND(rval).EQ.DOUBLE_KIND )then
1235	error = NF_PUT_ATT_DOUBLE( mpp_file(unit)%ncid, id, name, NF_BYTE, size(rval_i), rval )
1236	else if( KIND(rval).EQ.FLOAT_KIND )then
1237	error = NF_PUT_ATT_REAL ( mpp_file(unit)%ncid, id, name, NF_BYTE, size(rval_i), rval )
1238	end if
1239	call netcdf_err(error)
1240	deallocate(rval_i)
1241	else
1242	call mpp_error( FATAL, 'WRITE_ATTRIBUTE_NETCDF: only legal packing values are 1,2,4,8.' )
1243	end if
1244	else
1245	!default is to write FLOATs (32-bit)
1246	if( KIND(rval).EQ.DOUBLE_KIND )then
1247	error = NF_PUT_ATT_DOUBLE( mpp_file(unit)%ncid, id, name, NF_FLOAT, size(rval), rval )
1248	else if( KIND(rval).EQ.FLOAT_KIND )then
1249	error = NF_PUT_ATT_REAL ( mpp_file(unit)%ncid, id, name, NF_FLOAT, size(rval), rval )
1250	end if
1251	call netcdf_err(error)
1252	end if
1253	else if( PRESENT(ival) )then
1254	error = NF_PUT_ATT_INT ( mpp_file(unit)%ncid, id, name, NF_INT, size(ival), ival ); call netcdf_err(error)
1255	else if( present(cval) )then
1256	error = NF_PUT_ATT_TEXT( mpp_file(unit)%ncid, id, name, len_trim(cval), cval ); call netcdf_err(error)
1257	else
1258	call mpp_error( FATAL, 'WRITE_ATTRIBUTE_NETCDF: one of rval, ival, cval must be present.' )
1259	end if
1260	#endif /* use_netCDF */
1261	return
1262	end subroutine write_attribute_netcdf
1263
1264	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1265	! !
1266	! MPP_WRITE !
1267	! !
1268	! mpp_write is used to write data to the file on <unit> using the !
1269	! file parameters supplied by mpp_open(). Axis and field definitions !
1270	! must have previously been written to the file using mpp_write_meta. !
1271	! !
1272	! mpp_write can take 2 forms, one for distributed data and one for !
1273	! non-distributed data. Distributed data refer to arrays whose two !
1274	! fastest-varying indices are domain-decomposed. Distributed data !
1275	! must be 2D or 3D (in space). Non-distributed data can be 0-3D. !
1276	! !
1277	! In all calls to mpp_write, tstamp is an optional argument. It is to !
1278	! be omitted if the field was defined not to be a function of time. !
1279	! Results are unpredictable if the argument is supplied for a time- !
1280	! independent field, or omitted for a time-dependent field. Repeated !
1281	! writes of a time-independent field are also not recommended. One !
1282	! time level of one field is written per call. !
1283	! !
1284	! !
1285	! For non-distributed data, use !
1286	! !
1287	! mpp_write( unit, field, data, tstamp ) !
1288	! integer, intent(in) :: unit !
1289	! type(fieldtype), intent(in) :: field !
1290	! real, optional :: tstamp !
1291	! data is real and can be scalar or of rank 1-3. !
1292	! !
1293	! For distributed data, use !
1294	! !
1295	! mpp_write( unit, field, domain, data, tstamp ) !
1296	! integer, intent(in) :: unit !
1297	! type(fieldtype), intent(in) :: field !
1298	! type(domain2D), intent(in) :: domain !
1299	! real, optional :: tstamp !
1300	! data is real and can be of rank 2 or 3. !
1301	! !
1302	! mpp_write( unit, axis ) !
1303	! integer, intent(in) :: unit !
1304	! type(axistype), intent(in) :: axis !
1305	! !
1306	! This call writes the actual co-ordinate values along each space !
1307	! axis. It must be called once for each space axis after all other !
1308	! metadata has been written. !
1309	! !
1310	! The mpp_write package also includes the routine write_record which !
1311	! performs the actual write. This routine is private to this module. !
1312	! !
1313	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1314	#define MPP_WRITE_2DDECOMP_1D_ mpp_write_2ddecomp_r1d
1315	#define MPP_WRITE_2DDECOMP_2D_ mpp_write_2ddecomp_r2d
1316	#define MPP_WRITE_2DDECOMP_3D_ mpp_write_2ddecomp_r3d
1317	#define MPP_TYPE_ real
1318	#include <mpp_write_2Ddecomp.h>
1319
1320	#define MPP_WRITE_ mpp_write_r0D
1321	#define MPP_TYPE_ real
1322	#define MPP_RANK_ !
1323	#define MPP_WRITE_RECORD_ call write_record( unit, field, 1, (/data/), tstamp )
1324	#include <mpp_write.h>
1325
1326	#define MPP_WRITE_ mpp_write_r1D
1327	#define MPP_TYPE_ real
1328	#define MPP_WRITE_RECORD_ call write_record( unit, field, size(data), data, tstamp )
1329	#define MPP_RANK_ (:)
1330	#include <mpp_write.h>
1331
1332	#define MPP_WRITE_ mpp_write_r2D
1333	#define MPP_TYPE_ real
1334	#define MPP_WRITE_RECORD_ call write_record( unit, field, size(data), data, tstamp )
1335	#define MPP_RANK_ (:,:)
1336	#include <mpp_write.h>
1337
1338	#define MPP_WRITE_ mpp_write_r3D
1339	#define MPP_TYPE_ real
1340	#define MPP_WRITE_RECORD_ call write_record( unit, field, size(data), data, tstamp )
1341	#define MPP_RANK_ (:,:,:)
1342	#include <mpp_write.h>
1343
1344	#define MPP_WRITE_ mpp_write_r4D
1345	#define MPP_TYPE_ real
1346	#define MPP_WRITE_RECORD_ call write_record( unit, field, size(data), data, tstamp )
1347	#define MPP_RANK_ (:,:,:,:)
1348	#include <mpp_write.h>
1349
1350	subroutine mpp_write_axis( unit, axis )
1351	integer, intent(in) :: unit
1352	type(axistype), intent(in) :: axis
1353	type(fieldtype) :: field
1354	integer :: is, ie
1355
1356	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_WRITE: must first call mpp_io_init.' )
1357	if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_WRITE: invalid unit number.' )
1358	if( mpp_file(unit)%threading.EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
1359	if( mpp_file(unit)%fileset .EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
1360	!we convert axis to type(fieldtype) in order to call write_record
1361	field = default_field
1362	allocate( field%axes(1) )
1363	field%axes(1) = axis
1364	allocate( field%size(1) )
1365	field%id = axis%id
1366	if( mpp_file(unit)%fileset.EQ.MPP_MULTI .AND. axis%domain.NE.NULL_DOMAIN1D )then
1367	call mpp_get_compute_domain( axis%domain, is, ie )
1368	field%size(1) = ie-is+1
1369	call write_record( unit, field, field%size(1), axis%data(is:) )
1370	else
1371	field%size(1) = size(axis%data)
1372	call write_record( unit, field, field%size(1), axis%data )
1373	end if
1374	return
1375	end subroutine mpp_write_axis
1376
1377	subroutine write_record( unit, field, nwords, data, time_in, domain )
1378	!routine that is finally called by all mpp_write routines to perform the write
1379	!a non-netCDF record contains:
1380	! field ID
1381	! a set of 4 coordinates (is:ie,js:je) giving the data subdomain
1382	! a timelevel and a timestamp (=NULLTIME if field is static)
1383	! 3D real data (stored as 1D)
1384	!if you are using direct access I/O, the RECL argument to OPEN must be large enough for the above
1385	!in a global direct access file, record position on PE is given by %record.
1386
1387	!Treatment of timestamp:
1388	! We assume that static fields have been passed without a timestamp.
1389	! Here that is converted into a timestamp of NULLTIME.
1390	! For non-netCDF fields, field is treated no differently, but is written
1391	! with a timestamp of NULLTIME. There is no check in the code to prevent
1392	! the user from repeatedly writing a static field.
1393
1394	integer, intent(in) :: unit, nwords
1395	type(fieldtype), intent(in) :: field
1396	real, intent(in) :: data(nwords)
1397	real(DOUBLE_KIND), intent(in), optional :: time_in
1398	type(domain2D), intent(in), optional :: domain
1399	integer, dimension(size(field%axes)) :: start, axsiz
1400	real :: time
1401	integer :: time_level
1402	logical :: newtime
1403	integer :: subdomain(4)
1404	integer :: packed_data(nwords)
1405	integer :: i, is, ie, js, je, isg, ieg, jsg, jeg, isizc, jsizc, isizg, jsizg
1406
1407	#ifdef use_CRI_pointers
1408	real(FLOAT_KIND) :: data_r4(nwords)
1409	pointer( ptr1, data_r4)
1410	pointer( ptr2, packed_data)
1411
1412	if (mpp_io_stack_size < 2nwords) call mpp_io_set_stack_size(2nwords)
1413
1414	ptr1 = LOC(mpp_io_stack(1))
1415	ptr2 = LOC(mpp_io_stack(nwords+1))
1416	#endif
1417
1418	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_WRITE: must first call mpp_io_init.' )
1419	if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_WRITE: invalid unit number.' )
1420	if( mpp_file(unit)%threading.EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
1421	if( mpp_file(unit)%fileset .EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
1422
1423	if( .NOT.mpp_file(unit)%initialized )then
1424	!this is the first call to mpp_write
1425	!we now declare the file to be initialized
1426	!if this is netCDF we switch file from DEFINE mode to DATA mode
1427	if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
1428	#ifdef use_netCDF
1429	!NOFILL is probably required for parallel: any circumstances in which not advisable?
1430	error = NF_SET_FILL( mpp_file(unit)%ncid, NF_NOFILL, i ); call netcdf_err(error)
1431	if( mpp_file(unit)%action.EQ.MPP_WRONLY )error = NF_ENDDEF(mpp_file(unit)%ncid); call netcdf_err(error)
1432	#endif
1433	else
1434	call mpp_write_meta( unit, 'END', cval='metadata' )
1435	end if
1436	mpp_file(unit)%initialized = .TRUE.
1437	if( verbose )print '(a,i3,a)', 'MPP_WRITE: PE=', pe, ' initialized file '//trim(mpp_file(unit)%name)//'.'
1438	end if
1439
1440	!initialize time: by default assume NULLTIME
1441	time = NULLTIME
1442	time_level = -1
1443	newtime = .FALSE.
1444	if( PRESENT(time_in) )time = time_in
1445	!increment time level if new time
1446	if( time.GT.mpp_file(unit)%time+EPSILON(time) )then !new time
1447	mpp_file(unit)%time_level = mpp_file(unit)%time_level + 1
1448	mpp_file(unit)%time = time
1449	newtime = .TRUE.
1450	end if
1451	if( verbose )print '(a,2i3,2i5,es13.5)', 'MPP_WRITE: PE, unit, %id, %time_level, %time=',&
1452	pe, unit, mpp_file(unit)%id, mpp_file(unit)%time_level, mpp_file(unit)%time
1453
1454	if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
1455	!define netCDF data block to be written:
1456	! time axis: START = time level
1457	! AXSIZ = 1
1458	! space axis: if there is no domain info
1459	! START = 1
1460	! AXSIZ = field%size(axis)
1461	! if there IS domain info:
1462	! start of domain is compute%start_index for multi-file I/O
1463	! global%start_index for all other cases
1464	! this number must be converted to 1 for NF_PUT_VAR
1465	! (netCDF fortran calls are with reference to 1),
1466	! So, START = compute%start_index - <start of domain> + 1
1467	! AXSIZ = usually compute%size
1468	! However, if compute%start_index-compute%end_index+1.NE.compute%size,
1469	! we assume that the call is passing a subdomain.
1470	! To pass a subdomain, you must pass a domain2D object that satisfies the following:
1471	! global%start_index must contain the <start of domain> as defined above;
1472	! the data domain and compute domain must refer to the subdomain being passed.
1473	! In this case, START = compute%start_index - <start of domain> + 1
1474	! AXSIZ = compute%start_index - compute%end_index + 1
1475	! NOTE: passing of subdomains will fail for multi-PE single-threaded I/O,
1476	! since that attempts to gather all data on PE 0.
1477	start = 1
1478	do i = 1,size(field%axes)
1479	axsiz(i) = field%size(i)
1480	if( i.EQ.field%time_axis_index )start(i) = mpp_file(unit)%time_level
1481	start(i) = max(start(i),1)
1482	end do
1483	if( PRESENT(domain) )then
1484	call mpp_get_compute_domain( domain, is, ie, js, je, xsize=isizc, ysize=jsizc )
1485	call mpp_get_global_domain ( domain, isg, ieg, jsg, jeg, xsize=isizg, ysize=jsizg )
1486	axsiz(1) = isizc
1487	axsiz(2) = jsizc
1488	if( npes.GT.1 .AND. mpp_file(unit)%fileset.EQ.MPP_SINGLE )then
1489	start(1) = is - isg + 1
1490	start(2) = js - jsg + 1
1491	else
1492	if( isizc.NE.ie-is+1 )then
1493	start(1) = is - isg + 1
1494	axsiz(1) = ie - is + 1
1495	end if
1496	if( jsizc.NE.je-js+1 )then
1497	start(2) = js - jsg + 1
1498	axsiz(2) = je - js + 1
1499	end if
1500	end if
1501	end if
1502	if( debug )print '(a,2i3,12i4)', 'WRITE_RECORD: PE, unit, start, axsiz=', pe, unit, start, axsiz
1503	#ifdef use_netCDF
1504	!write time information if new time
1505	if( newtime )then
1506	if( KIND(time).EQ.DOUBLE_KIND )then
1507	error = NF_PUT_VAR1_DOUBLE( mpp_file(unit)%ncid, mpp_file(unit)%id, mpp_file(unit)%time_level, time )
1508	else if( KIND(time).EQ.FLOAT_KIND )then
1509	error = NF_PUT_VAR1_REAL ( mpp_file(unit)%ncid, mpp_file(unit)%id, mpp_file(unit)%time_level, time )
1510	end if
1511	end if
1512	if( field%pack.LE.2 )then
1513	if( KIND(data).EQ.DOUBLE_KIND )then
1514	! write(stderr,*)data
1515	error = NF_PUT_VARA_DOUBLE( mpp_file(unit)%ncid, field%id, start, axsiz, data )
1516	else if( KIND(data).EQ.FLOAT_KIND )then
1517	error = NF_PUT_VARA_REAL ( mpp_file(unit)%ncid, field%id, start, axsiz, data )
1518	end if
1519	else !convert to integer using scale and add: no error check on packed data representation
1520	packed_data = nint((data-field%add)/field%scale)
1521	error = NF_PUT_VARA_INT ( mpp_file(unit)%ncid, field%id, start, axsiz, packed_data )
1522	end if
1523	call netcdf_err(error)
1524	#endif
1525	else !non-netCDF
1526	!subdomain contains (/is,ie,js,je/)
1527	if( PRESENT(domain) )then
1528	subdomain(:) = (/ is, ie, js, je /)
1529	else
1530	subdomain(:) = -1 ! -1 means use global value from axis metadata
1531	end if
1532	if( mpp_file(unit)%format.EQ.MPP_ASCII )then
1533	!implies sequential access
1534	write( unit,* )field%id, subdomain, time_level, time, data
1535	else !MPP_IEEE32 or MPP_NATIVE
1536	if( mpp_file(unit)%access.EQ.MPP_SEQUENTIAL )then
1537	#ifdef __sgi
1538	if( mpp_file(unit)%format.EQ.MPP_IEEE32 )then
1539	data_r4 = data !IEEE conversion layer on SGI until assign -N ieee_32 is supported
1540	write(unit)field%id, subdomain, time_level, time, data_r4
1541	else
1542	write(unit)field%id, subdomain, time_level, time, data
1543	end if
1544	#else
1545	write(unit)field%id, subdomain, time_level, time, data
1546	#endif
1547	else !MPP_DIRECT
1548	#ifdef __sgi
1549	if( mpp_file(unit)%format.EQ.MPP_IEEE32 )then
1550	data_r4 = data !IEEE conversion layer on SGI until assign -N ieee_32 is supported
1551	write( unit, rec=mpp_file(unit)%record )field%id, subdomain, time_level, time, data_r4
1552	else
1553	write( unit, rec=mpp_file(unit)%record )field%id, subdomain, time_level, time, data
1554	end if
1555	#else
1556	write( unit, rec=mpp_file(unit)%record )field%id, subdomain, time_level, time, data
1557	#endif
1558	if( debug )print '(a,i3,a,i3)', 'MPP_WRITE: PE=', pe, ' wrote record ', mpp_file(unit)%record
1559	end if
1560	end if
1561	end if
1562
1563	!recompute current record for direct access I/O
1564	if( mpp_file(unit)%access.EQ.MPP_DIRECT )then
1565	if( mpp_file(unit)%fileset.EQ.MPP_SINGLE )then
1566	!assumes all PEs participate in I/O: modify later
1567	mpp_file(unit)%record = mpp_file(unit)%record + records_per_pe*npes
1568	else
1569	mpp_file(unit)%record = mpp_file(unit)%record + records_per_pe
1570	end if
1571	end if
1572
1573	return
1574	end subroutine write_record
1575
1576	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1577	! !
1578	! MPP_COPY_META !
1579	! !
1580	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1581	subroutine mpp_copy_meta_global( unit, gatt )
1582	!writes a global metadata attribute to unit <unit>
1583	!attribute <name> can be an real, integer or character
1584	!one and only one of rval, ival, and cval should be present
1585	!the first found will be used
1586	!for a non-netCDF file, it is encoded into a string "GLOBAL <name> <val>"
1587	integer, intent(in) :: unit
1588	type(atttype), intent(in) :: gatt
1589	integer :: len
1590
1591	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_WRITE_META: must first call mpp_io_init.' )
1592	if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_WRITE_META: invalid unit number.' )
1593	if( mpp_file(unit)%threading.EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
1594	if( mpp_file(unit)%fileset.EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
1595	if( mpp_file(unit)%action.NE.MPP_WRONLY )return !no writing metadata on APPEND
1596	if( mpp_file(unit)%initialized ) &
1597	call mpp_error( FATAL, 'MPP_WRITE_META: cannot write metadata to file after an mpp_write.' )
1598	#ifdef use_netCDF
1599	if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
1600	if( gatt%type.EQ.NF_CHAR )then
1601	len = gatt%len
1602	call write_attribute_netcdf( unit, NF_GLOBAL, gatt%name, cval=gatt%catt(1:len) )
1603	else
1604	call write_attribute_netcdf( unit, NF_GLOBAL, gatt%name, rval=gatt%fatt )
1605	endif
1606	else
1607	if( gatt%type.EQ.NF_CHAR )then
1608	len=gatt%len
1609	call write_attribute( unit, 'GLOBAL '//trim(gatt%name), cval=gatt%catt(1:len) )
1610	else
1611	call write_attribute( unit, 'GLOBAL '//trim(gatt%name), rval=gatt%fatt )
1612	endif
1613	end if
1614	#else
1615	call mpp_error( FATAL, 'MPP_READ currently requires use_netCDF option' )
1616	#endif
1617	return
1618	end subroutine mpp_copy_meta_global
1619
1620	subroutine mpp_copy_meta_axis( unit, axis, domain )
1621	!load the values in an axistype (still need to call mpp_write)
1622	!write metadata attributes for axis. axis is declared inout
1623	!because the variable and dimension ids are altered
1624
1625	integer, intent(in) :: unit
1626	type(axistype), intent(inout) :: axis
1627	type(domain1D), intent(in), optional :: domain
1628	character(len=512) :: text
1629	integer :: i, len, is, ie, isg, ieg
1630
1631	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_WRITE_META: must first call mpp_io_init.' )
1632	if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_WRITE_META: invalid unit number.' )
1633	if( mpp_file(unit)%threading.EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
1634	if( mpp_file(unit)%fileset.EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
1635	if( mpp_file(unit)%action.NE.MPP_WRONLY )return !no writing metadata on APPEND
1636	if( mpp_file(unit)%initialized ) &
1637	call mpp_error( FATAL, 'MPP_WRITE_META: cannot write metadata to file after an mpp_write.' )
1638
1639	! redefine domain if present
1640	if( PRESENT(domain) )then
1641	axis%domain = domain
1642	else
1643	axis%domain = NULL_DOMAIN1D
1644	end if
1645
1646	#ifdef use_netCDF
1647	!write metadata
1648	if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
1649
1650	!write axis def
1651	if( ASSOCIATED(axis%data) )then !space axis
1652	if( mpp_file(unit)%fileset.EQ.MPP_MULTI .AND. axis%domain.NE.NULL_DOMAIN1D )then
1653	call mpp_get_compute_domain( axis%domain, is, ie )
1654	call mpp_get_global_domain( axis%domain, isg, ieg )
1655	error = NF_DEF_DIM( mpp_file(unit)%ncid, axis%name, ie-is+1, axis%did )
1656	else
1657	error = NF_DEF_DIM( mpp_file(unit)%ncid, axis%name, size(axis%data), axis%did )
1658	end if
1659	call netcdf_err(error)
1660	error = NF_DEF_VAR( mpp_file(unit)%ncid, axis%name, NF_FLOAT, 1, axis%did, axis%id ); call netcdf_err(error)
1661	else !time axis
1662	error = NF_DEF_DIM( mpp_file(unit)%ncid, axis%name, NF_UNLIMITED, axis%did ); call netcdf_err(error)
1663	error = NF_DEF_VAR( mpp_file(unit)%ncid, axis%name, NF_DOUBLE, 1, axis%did, axis%id ); call netcdf_err(error)
1664	mpp_file(unit)%id = axis%id !file ID is the same as time axis varID
1665	mpp_file(unit)%recdimid = axis%did ! record dimension id
1666	end if
1667	else
1668	varnum = varnum + 1
1669	axis%id = varnum
1670	axis%did = varnum
1671	!write axis def
1672	write( text, '(a,i4,a)' )'AXIS ', axis%id, ' name'
1673	call write_attribute( unit, trim(text), cval=axis%name )
1674	write( text, '(a,i4,a)' )'AXIS ', axis%id, ' size'
1675	if( ASSOCIATED(axis%data) )then !space axis
1676	if( mpp_file(unit)%fileset.EQ.MPP_MULTI .AND. axis%domain.NE.NULL_DOMAIN1D )then
1677	call write_attribute( unit, trim(text), ival=(/ie-is+1/) )
1678	else
1679	call write_attribute( unit, trim(text), ival=(/size(axis%data)/) )
1680	end if
1681	else !time axis
1682	if( mpp_file(unit)%id.NE.-1 ) &
1683	call mpp_error( FATAL, 'MPP_WRITE_META_AXIS: There is already a time axis for this file.' )
1684	call write_attribute( unit, trim(text), ival=(/0/) ) !a size of 0 indicates time axis
1685	mpp_file(unit)%id = axis%id
1686	end if
1687	end if
1688	!write axis attributes
1689
1690	do i=1,axis%natt
1691	if( axis%Att(i)%name.NE.default_att%name )then
1692	if( axis%Att(i)%type.EQ.NF_CHAR )then
1693	len = axis%Att(i)%len
1694	call mpp_write_meta( unit, axis%id, axis%Att(i)%name, cval=axis%Att(i)%catt(1:len) )
1695	else
1696	call mpp_write_meta( unit, axis%id, axis%Att(i)%name, rval=axis%Att(i)%fatt)
1697	endif
1698	endif
1699	enddo
1700
1701	if( mpp_file(unit)%threading.EQ.MPP_MULTI .AND. mpp_file(unit)%fileset.EQ.MPP_MULTI .AND. axis%domain.NE.NULL_DOMAIN1D )then
1702	call mpp_write_meta( unit, axis%id, 'domain_decomposition', ival=(/isg,ieg,is,ie/) )
1703	end if
1704	if( verbose )print '(a,2i3,x,a,2i3)', 'MPP_WRITE_META: Wrote axis metadata, pe, unit, axis%name, axis%id, axis%did=', &
1705	pe, unit, trim(axis%name), axis%id, axis%did
1706	#else
1707	call mpp_error( FATAL, 'MPP_READ currently requires use_netCDF option' )
1708	#endif
1709	return
1710	end subroutine mpp_copy_meta_axis
1711
1712	subroutine mpp_copy_meta_field( unit, field, axes )
1713	!useful for copying field metadata from a previous call to mpp_read_meta
1714	!define field: must have already called mpp_write_meta(axis) for each axis
1715	integer, intent(in) :: unit
1716	type(fieldtype), intent(inout) :: field
1717	type(axistype), intent(in), optional :: axes(:)
1718	!this array is required because of f77 binding on netCDF interface
1719	integer, allocatable :: axis_id(:)
1720	real :: a, b
1721	integer :: i
1722
1723	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_WRITE_META: must first call mpp_io_init.' )
1724	if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_WRITE_META: invalid unit number.' )
1725	if( mpp_file(unit)%threading.EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
1726	if( mpp_file(unit)%fileset.EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
1727	if( mpp_file(unit)%action.NE.MPP_WRONLY )return !no writing metadata on APPEND
1728	if( mpp_file(unit)%initialized ) &
1729	call mpp_error( FATAL, 'MPP_WRITE_META: cannot write metadata to file after an mpp_write.' )
1730
1731	if( field%pack.NE.1 .AND. field%pack.NE.2 )then
1732	if( field%pack.NE.4 .AND. field%pack.NE.8 ) &
1733	call mpp_error( FATAL, 'MPP_WRITE_META_FIELD: only legal packing values are 1,2,4,8.' )
1734	end if
1735
1736	if (PRESENT(axes)) then
1737	deallocate(field%axes)
1738	deallocate(field%size)
1739	allocate(field%axes(size(axes)))
1740	allocate(field%size(size(axes)))
1741	field%axes = axes
1742	do i=1,size(axes)
1743	if (ASSOCIATED(axes(i)%data)) then
1744	field%size(i) = size(axes(i)%data)
1745	else
1746	field%size(i) = 1
1747	field%time_axis_index = i
1748	endif
1749	enddo
1750	endif
1751
1752	if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
1753	#ifdef use_netCDF
1754	allocate( axis_id(size(field%axes)) )
1755	do i = 1,size(field%axes)
1756	axis_id(i) = field%axes(i)%did
1757	end do
1758	!write field def
1759	select case (field%pack)
1760	case(1)
1761	error = NF_DEF_VAR( mpp_file(unit)%ncid, field%name, NF_DOUBLE, size(field%axes), axis_id, field%id )
1762	case(2)
1763	error = NF_DEF_VAR( mpp_file(unit)%ncid, field%name, NF_FLOAT, size(field%axes), axis_id, field%id )
1764	case(4)
1765	if( field%scale.EQ.default_field%scale .OR. field%add.EQ.default_field%add ) &
1766	call mpp_error( FATAL, 'MPP_WRITE_META_FIELD: scale and add must be supplied when pack=4.' )
1767	error = NF_DEF_VAR( mpp_file(unit)%ncid, field%name, NF_SHORT, size(field%axes), axis_id, field%id )
1768	case(8)
1769	if( field%scale.EQ.default_field%scale .OR. field%add.EQ.default_field%add ) &
1770	call mpp_error( FATAL, 'MPP_WRITE_META_FIELD: scale and add must be supplied when pack=8.' )
1771	error = NF_DEF_VAR( mpp_file(unit)%ncid, field%name, NF_BYTE, size(field%axes), axis_id, field%id )
1772	case default
1773	call mpp_error( FATAL, 'MPP_WRITE_META_FIELD: only legal packing values are 1,2,4,8.' )
1774	end select
1775	#endif
1776	else
1777	varnum = varnum + 1
1778	field%id = varnum
1779	if( field%pack.NE.default_field%pack ) &
1780	call mpp_error( WARNING, 'MPP_WRITE_META: Packing is currently available only on netCDF files.' )
1781	!write field def
1782	write( text, '(a,i4,a)' )'FIELD ', field%id, ' name'
1783	call write_attribute( unit, trim(text), cval=field%name )
1784	write( text, '(a,i4,a)' )'FIELD ', field%id, ' axes'
1785	call write_attribute( unit, trim(text), ival=field%axes(:)%did )
1786	end if
1787	!write field attributes: these names follow netCDF conventions
1788	call mpp_write_meta( unit, field%id, 'long_name', cval=field%longname )
1789	call mpp_write_meta( unit, field%id, 'units', cval=field%units )
1790	!all real attributes must be written as packed
1791	if( (field%min.NE.default_field%min) .AND. (field%max.NE.default_field%max) )then
1792	if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
1793	call mpp_write_meta( unit, field%id, 'valid_range', rval=(/field%min,field%max/), pack=field%pack )
1794	else
1795	a = nint((field%min-field%add)/field%scale)
1796	b = nint((field%max-field%add)/field%scale)
1797	call mpp_write_meta( unit, field%id, 'valid_range', rval=(/a, b /), pack=field%pack )
1798	end if
1799	else if( field%min.NE.default_field%min )then
1800	if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
1801	call mpp_write_meta( unit, field%id, 'valid_min', rval=field%min, pack=field%pack )
1802	else
1803	a = nint((field%min-field%add)/field%scale)
1804	call mpp_write_meta( unit, field%id, 'valid_min', rval=a, pack=field%pack )
1805	end if
1806	else if( field%max.NE.default_field%max )then
1807	if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
1808	call mpp_write_meta( unit, field%id, 'valid_max', rval=field%max, pack=field%pack )
1809	else
1810	a = nint((field%max-field%add)/field%scale)
1811	call mpp_write_meta( unit, field%id, 'valid_max', rval=a, pack=field%pack )
1812	end if
1813	end if
1814	if( field%missing.NE.default_field%missing )then
1815	if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
1816	call mpp_write_meta( unit, field%id, 'missing_value', rval=field%missing, pack=field%pack )
1817	else
1818	a = nint((field%missing-field%add)/field%scale)
1819	call mpp_write_meta( unit, field%id, 'missing_value', rval=a, pack=field%pack )
1820	end if
1821	end if
1822	if( field%fill.NE.default_field%fill )then
1823	if( field%pack.EQ.1 .OR. field%pack.EQ.2 )then
1824	call mpp_write_meta( unit, field%id, '_FillValue', rval=field%missing, pack=field%pack )
1825	else
1826	a = nint((field%fill-field%add)/field%scale)
1827	call mpp_write_meta( unit, field%id, '_FillValue', rval=a, pack=field%pack )
1828	end if
1829	end if
1830	if( field%pack.NE.1 .AND. field%pack.NE.2 )then
1831	call mpp_write_meta( unit, field%id, 'packing', ival=field%pack )
1832	if( field%scale.NE.default_field%scale )call mpp_write_meta( unit, field%id, 'scale_factor', rval=field%scale )
1833	if( field%add.NE.default_field%add )call mpp_write_meta( unit, field%id, 'add_offset', rval=field%add )
1834	end if
1835	if( verbose )print '(a,2i3,x,a,i3)', 'MPP_WRITE_META: Wrote field metadata: pe, unit, field%name, field%id=', &
1836	pe, unit, trim(field%name), field%id
1837
1838	return
1839	end subroutine mpp_copy_meta_field
1840
1841	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1842	! !
1843	! MPP_READ !
1844	! !
1845	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1846
1847	#define MPP_READ_2DDECOMP_1D_ mpp_read_2ddecomp_r1d
1848	#define MPP_READ_2DDECOMP_2D_ mpp_read_2ddecomp_r2d
1849	#define MPP_READ_2DDECOMP_3D_ mpp_read_2ddecomp_r3d
1850	#define MPP_TYPE_ real
1851	#include <mpp_read_2Ddecomp.h>
1852
1853	subroutine read_record( unit, field, nwords, data, time_level, domain )
1854	!routine that is finally called by all mpp_read routines to perform the read
1855	!a non-netCDF record contains:
1856	! field ID
1857	! a set of 4 coordinates (is:ie,js:je) giving the data subdomain
1858	! a timelevel and a timestamp (=NULLTIME if field is static)
1859	! 3D real data (stored as 1D)
1860	!if you are using direct access I/O, the RECL argument to OPEN must be large enough for the above
1861	!in a global direct access file, record position on PE is given by %record.
1862
1863	!Treatment of timestamp:
1864	! We assume that static fields have been passed without a timestamp.
1865	! Here that is converted into a timestamp of NULLTIME.
1866	! For non-netCDF fields, field is treated no differently, but is written
1867	! with a timestamp of NULLTIME. There is no check in the code to prevent
1868	! the user from repeatedly writing a static field.
1869
1870	integer, intent(in) :: unit, nwords
1871	type(fieldtype), intent(in) :: field
1872	real, intent(inout) :: data(nwords)
1873	integer, intent(in), optional :: time_level
1874	type(domain2D), intent(in), optional :: domain
1875	integer, dimension(size(field%axes)) :: start, axsiz
1876	real :: time
1877
1878	logical :: newtime
1879	integer :: subdomain(4), tlevel
1880
1881	integer(SHORT_KIND) :: i2vals(nwords)
1882	!#ifdef __sgi
1883	integer(INT_KIND) :: ivals(nwords)
1884	real(FLOAT_KIND) :: rvals(nwords)
1885	!#else
1886	! integer :: ivals(nwords)
1887	! real :: rvals(nwords)
1888	!#endif
1889
1890	real(DOUBLE_KIND) :: r8vals(nwords)
1891
1892	integer :: i, error, is, ie, js, je, isg, ieg, jsg, jeg
1893
1894	#ifdef use_CRI_pointers
1895	pointer( ptr1, i2vals )
1896	pointer( ptr2, ivals )
1897	pointer( ptr3, rvals )
1898	pointer( ptr4, r8vals )
1899
1900	if (mpp_io_stack_size < 4nwords) call mpp_io_set_stack_size(4nwords)
1901
1902	ptr1 = LOC(mpp_io_stack(1))
1903	ptr2 = LOC(mpp_io_stack(nwords+1))
1904	ptr3 = LOC(mpp_io_stack(2*nwords+1))
1905	ptr4 = LOC(mpp_io_stack(3*nwords+1))
1906	#endif
1907	if (.not.PRESENT(time_level)) then
1908	tlevel = 0
1909	else
1910	tlevel = time_level
1911	endif
1912
1913	#ifdef use_netCDF
1914	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'READ_RECORD: must first call mpp_io_init.' )
1915	if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'READ_RECORD: invalid unit number.' )
1916	if( mpp_file(unit)%threading.EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
1917	if( mpp_file(unit)%fileset.EQ.MPP_MULTI )call mpp_error( FATAL, 'READ_RECORD: multiple filesets not supported for MPP_READ' )
1918
1919	if( .NOT.mpp_file(unit)%initialized ) call mpp_error( FATAL, 'MPP_READ: must first call mpp_read_meta.' )
1920
1921
1922
1923	if( verbose )print '(a,2i3,2i5)', 'MPP_READ: PE, unit, %id, %time_level =',&
1924	pe, unit, mpp_file(unit)%id, tlevel
1925
1926	if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
1927	!define netCDF data block to be read:
1928	! time axis: START = time level
1929	! AXSIZ = 1
1930	! space axis: if there is no domain info
1931	! START = 1
1932	! AXSIZ = field%size(axis)
1933	! if there IS domain info:
1934	! start of domain is compute%start_index for multi-file I/O
1935	! global%start_index for all other cases
1936	! this number must be converted to 1 for NF_GET_VAR
1937	! (netCDF fortran calls are with reference to 1),
1938	! So, START = compute%start_index - <start of domain> + 1
1939	! AXSIZ = usually compute%size
1940	! However, if compute%start_index-compute%end_index+1.NE.compute%size,
1941	! we assume that the call is passing a subdomain.
1942	! To pass a subdomain, you must pass a domain2D object that satisfies the following:
1943	! global%start_index must contain the <start of domain> as defined above;
1944	! the data domain and compute domain must refer to the subdomain being passed.
1945	! In this case, START = compute%start_index - <start of domain> + 1
1946	! AXSIZ = compute%start_index - compute%end_index + 1
1947	! NOTE: passing of subdomains will fail for multi-PE single-threaded I/O,
1948	! since that attempts to gather all data on PE 0.
1949	start = 1
1950	do i = 1,size(field%axes)
1951	axsiz(i) = field%size(i)
1952	if( field%axes(i)%did.EQ.field%time_axis_index )start(i) = tlevel
1953	end do
1954	if( PRESENT(domain) )then
1955	call mpp_get_compute_domain( domain, is, ie, js, je )
1956	call mpp_get_global_domain ( domain, isg, ieg, jsg, jeg )
1957	axsiz(1) = ie-is+1
1958	axsiz(2) = je-js+1
1959	if( npes.GT.1 .AND. mpp_file(unit)%fileset.EQ.MPP_SINGLE )then
1960	start(1) = is - isg + 1
1961	start(2) = js - jsg + 1
1962	else
1963	if( ie-is+1.NE.ie-is+1 )then
1964	start(1) = is - isg + 1
1965	axsiz(1) = ie - is + 1
1966	end if
1967	if( je-js+1.NE.je-js+1 )then
1968	start(2) = js - jsg + 1
1969	axsiz(2) = je - js + 1
1970	end if
1971	end if
1972	end if
1973
1974	if( verbose )print '(a,2i3,i6,12i4)', 'READ_RECORD: PE, unit, nwords, start, axsiz=', pe, unit, nwords, start, axsiz
1975
1976	select case (field%type)
1977	case(NF_BYTE)
1978	! use type conversion
1979	call mpp_error( FATAL, 'MPP_READ: does not support NF_BYTE packing' )
1980	case(NF_SHORT)
1981	error = NF_GET_VARA_INT2 ( mpp_file(unit)%ncid, field%id, start, axsiz, i2vals ); call netcdf_err(error)
1982	data(:)=i2vals(:)*field%scale + field%add
1983	case(NF_INT)
1984	error = NF_GET_VARA_INT ( mpp_file(unit)%ncid, field%id, start, axsiz, ivals ); call netcdf_err(error)
1985	data(:)=ivals(:)
1986	case(NF_FLOAT)
1987	error = NF_GET_VARA_REAL ( mpp_file(unit)%ncid, field%id, start, axsiz, rvals ); call netcdf_err(error)
1988	data(:)=rvals(:)
1989	case(NF_DOUBLE)
1990	error = NF_GET_VARA_DOUBLE( mpp_file(unit)%ncid, field%id, start, axsiz, r8vals ); call netcdf_err(error)
1991	data(:)=r8vals(:)
1992	case default
1993	call mpp_error( FATAL, 'MPP_READ: invalid pack value' )
1994	end select
1995	else !non-netCDF
1996	!subdomain contains (/is,ie,js,je/)
1997	call mpp_error( FATAL, 'Currently dont support non-NetCDF mpp read' )
1998
1999	end if
2000	#else
2001	call mpp_error( FATAL, 'MPP_READ currently requires use_netCDF option' )
2002	#endif
2003	return
2004	end subroutine read_record
2005
2006	subroutine mpp_read_r3D( unit, field, data, tindex)
2007	integer, intent(in) :: unit
2008	type(fieldtype), intent(in) :: field
2009	real, intent(inout) :: data(:,:,:)
2010	integer, intent(in), optional :: tindex
2011
2012	call read_record( unit, field, size(data), data, tindex )
2013	end subroutine mpp_read_r3D
2014
2015	subroutine mpp_read_r2D( unit, field, data, tindex )
2016	integer, intent(in) :: unit
2017	type(fieldtype), intent(in) :: field
2018	real, intent(inout) :: data(:,:)
2019	integer, intent(in), optional :: tindex
2020
2021	call read_record( unit, field, size(data), data, tindex )
2022	end subroutine mpp_read_r2D
2023
2024	subroutine mpp_read_r1D( unit, field, data, tindex )
2025	integer, intent(in) :: unit
2026	type(fieldtype), intent(in) :: field
2027	real, intent(inout) :: data(:)
2028	integer, intent(in), optional :: tindex
2029
2030	call read_record( unit, field, size(data), data, tindex )
2031	end subroutine mpp_read_r1D
2032
2033	subroutine mpp_read_r0D( unit, field, data, tindex )
2034	integer, intent(in) :: unit
2035	type(fieldtype), intent(in) :: field
2036	real, intent(inout) :: data
2037	integer, intent(in), optional :: tindex
2038	real, dimension(1) :: data_tmp
2039
2040	data_tmp(1)=data
2041	call read_record( unit, field, 1, data_tmp, tindex )
2042	data=data_tmp(1)
2043	end subroutine mpp_read_r0D
2044
2045	subroutine mpp_read_meta(unit)
2046	!
2047	! read file attributes including dimension and variable attributes
2048	! and store in filetype structure. All of the file information
2049	! with the exception of the (variable) data is stored. Attributes
2050	! are supplied to the user by get_info,get_atts,get_axes and get_fields
2051	!
2052	! every PE is eligible to call mpp_read_meta
2053	!
2054	integer, parameter :: MAX_DIMVALS = 100000
2055	integer, intent(in) :: unit
2056
2057	integer :: ncid,ndim,nvar_total,natt,recdim,nv,nvar,len
2058	integer :: error,i,j
2059	integer :: type,nvdims,nvatts, dimid
2060	integer, allocatable, dimension(:) :: dimids
2061	type(axistype) , allocatable, dimension(:) :: Axis
2062	character(len=128) :: name, attname, unlimname, attval
2063	logical :: isdim
2064
2065	integer(SHORT_KIND) :: i2vals(MAX_DIMVALS)
2066	!#ifdef __sgi
2067	integer(INT_KIND) :: ivals(MAX_DIMVALS)
2068	real(FLOAT_KIND) :: rvals(MAX_DIMVALS)
2069	!#else
2070	! integer :: ivals(MAX_DIMVALS)
2071	! real :: rvals(MAX_DIMVALS)
2072	!#endif
2073	real(DOUBLE_KIND) :: r8vals(MAX_DIMVALS)
2074
2075	#ifdef use_netCDF
2076
2077	if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
2078	ncid = mpp_file(unit)%ncid
2079	error = NF_INQ(ncid,ndim, nvar_total,&
2080	natt, recdim);call netcdf_err(error)
2081
2082
2083	mpp_file(unit)%ndim = ndim
2084	mpp_file(unit)%natt = natt
2085	mpp_file(unit)%recdimid = recdim
2086	!
2087	! if no recdim exists, recdimid = -1
2088	! variable id of unlimdim and length
2089	!
2090	if( recdim.NE.-1 )then
2091	error = NF_INQ_DIM( ncid, recdim, unlimname, mpp_file(unit)%time_level );call netcdf_err(error)
2092	error = NF_INQ_VARID( ncid, unlimname, mpp_file(unit)%id ); call netcdf_err(error)
2093	else
2094	mpp_file(unit)%time_level = -1 ! set to zero so mpp_get_info returns ntime=0 if no time axis present
2095	endif
2096
2097	allocate(mpp_file(unit)%Att(natt))
2098	allocate(Axis(ndim))
2099	allocate(dimids(ndim))
2100	allocate(mpp_file(unit)%Axis(ndim))
2101
2102	!
2103	! initialize fieldtype and axis type
2104	!
2105
2106
2107	do i=1,ndim
2108	Axis(i) = default_axis
2109	mpp_file(unit)%Axis(i) = default_axis
2110	enddo
2111
2112	do i=1,natt
2113	mpp_file(unit)%Att(i) = default_att
2114	enddo
2115
2116	!
2117	! assign global attributes
2118	!
2119	do i=1,natt
2120	error=NF_INQ_ATTNAME(ncid,NF_GLOBAL,i,name);call netcdf_err(error)
2121	error=NF_INQ_ATT(ncid,NF_GLOBAL,trim(name),type,len);call netcdf_err(error)
2122	mpp_file(unit)%Att(i)%name = name
2123	mpp_file(unit)%Att(i)%len = len
2124	mpp_file(unit)%Att(i)%type = type
2125	!
2126	! allocate space for att data and assign
2127	!
2128	select case (type)
2129	case (NF_CHAR)
2130	if (len.gt.512) then
2131	call mpp_error(NOTE,'GLOBAL ATT too long - not reading this metadata')
2132	len=7
2133	mpp_file(unit)%Att(i)%len=len
2134	mpp_file(unit)%Att(i)%catt = 'unknown'
2135	else
2136	error=NF_GET_ATT_TEXT(ncid,NF_GLOBAL,name,mpp_file(unit)%Att(i)%catt);call netcdf_err(error)
2137	if (verbose.and.pe == 0) print *, 'GLOBAL ATT ',trim(name),' ',mpp_file(unit)%Att(i)%catt(1:len)
2138	endif
2139	!
2140	! store integers in float arrays
2141	!
2142	case (NF_SHORT)
2143	allocate(mpp_file(unit)%Att(i)%fatt(len))
2144	error=NF_GET_ATT_INT2(ncid,NF_GLOBAL,name,i2vals);call netcdf_err(error)
2145	if( verbose .and. pe == 0 )print *, 'GLOBAL ATT ',trim(name),' ',i2vals(1:len)
2146	mpp_file(unit)%Att(i)%fatt(1:len)=i2vals(1:len)
2147	case (NF_INT)
2148	allocate(mpp_file(unit)%Att(i)%fatt(len))
2149	error=NF_GET_ATT_INT(ncid,NF_GLOBAL,name,ivals);call netcdf_err(error)
2150	if( verbose .and. pe == 0 )print *, 'GLOBAL ATT ',trim(name),' ',ivals(1:len)
2151	mpp_file(unit)%Att(i)%fatt(1:len)=ivals(1:len)
2152	case (NF_FLOAT)
2153	allocate(mpp_file(unit)%Att(i)%fatt(len))
2154	error=NF_GET_ATT_REAL(ncid,NF_GLOBAL,name,rvals);call netcdf_err(error)
2155	mpp_file(unit)%Att(i)%fatt(1:len)=rvals(1:len)
2156	if( verbose .and. pe == 0)print *, 'GLOBAL ATT ',trim(name),' ',mpp_file(unit)%Att(i)%fatt(1:len)
2157	case (NF_DOUBLE)
2158	allocate(mpp_file(unit)%Att(i)%fatt(len))
2159	error=NF_GET_ATT_DOUBLE(ncid,NF_GLOBAL,name,r8vals);call netcdf_err(error)
2160	mpp_file(unit)%Att(i)%fatt(1:len)=r8vals(1:len)
2161	if( verbose .and. pe == 0)print *, 'GLOBAL ATT ',trim(name),' ',mpp_file(unit)%Att(i)%fatt(1:len)
2162	end select
2163
2164	enddo
2165	!
2166	! assign dimension name and length
2167	!
2168	do i=1,ndim
2169	error = NF_INQ_DIM(ncid,i,name,len);call netcdf_err(error)
2170	Axis(i)%name = name
2171	Axis(i)%len = len
2172	enddo
2173
2174	nvar=0
2175	do i=1, nvar_total
2176	error=NF_INQ_VAR(ncid,i,name,type,nvdims,dimids,nvatts);call netcdf_err(error)
2177	isdim=.false.
2178	do j=1,ndim
2179	if( trim(lowercase(name)).EQ.trim(lowercase(Axis(j)%name)) )isdim=.true.
2180	enddo
2181	if (.not.isdim) nvar=nvar+1
2182	enddo
2183	mpp_file(unit)%nvar = nvar
2184	allocate(mpp_file(unit)%Var(nvar))
2185
2186	do i=1,nvar
2187	mpp_file(unit)%Var(i) = default_field
2188	enddo
2189
2190	!
2191	! assign dimension info
2192	!
2193	do i=1, nvar_total
2194	error=NF_INQ_VAR(ncid,i,name,type,nvdims,dimids,nvatts);call netcdf_err(error)
2195	isdim=.false.
2196	do j=1,ndim
2197	if( trim(lowercase(name)).EQ.trim(lowercase(Axis(j)%name)) )isdim=.true.
2198	enddo
2199
2200	if( isdim )then
2201	error=NF_INQ_DIMID(ncid,name,dimid);call netcdf_err(error)
2202	Axis(dimid)%type = type
2203	Axis(dimid)%did = dimid
2204	Axis(dimid)%id = i
2205	Axis(dimid)%natt = nvatts
2206	! get axis values
2207	if( i.NE.mpp_file(unit)%id )then ! non-record dims
2208	select case (type)
2209	case (NF_INT)
2210	len=Axis(dimid)%len
2211	allocate(Axis(dimid)%data(len))
2212	error = NF_GET_VAR_INT(ncid,i,ivals);call netcdf_err(error)
2213	Axis(dimid)%data(1:len)=ivals(1:len)
2214	case (NF_FLOAT)
2215	len=Axis(dimid)%len
2216	allocate(Axis(dimid)%data(len))
2217	error = NF_GET_VAR_REAL(ncid,i,rvals);call netcdf_err(error)
2218	Axis(dimid)%data(1:len)=rvals(1:len)
2219	case (NF_DOUBLE)
2220	len=Axis(dimid)%len
2221	allocate(Axis(dimid)%data(len))
2222	error = NF_GET_VAR_DOUBLE(ncid,i,r8vals);call netcdf_err(error)
2223	Axis(dimid)%data(1:len) = r8vals(1:len)
2224	case default
2225	call mpp_error( FATAL, 'Invalid data type for dimension' )
2226	end select
2227	else
2228	len = mpp_file(unit)%time_level
2229	allocate(mpp_file(unit)%time_values(len))
2230	select case (type)
2231	case (NF_FLOAT)
2232	error = NF_GET_VAR_REAL(ncid,i,rvals);call netcdf_err(error)
2233	mpp_file(unit)%time_values(1:len) = rvals(1:len)
2234	case (NF_DOUBLE)
2235	error = NF_GET_VAR_DOUBLE(ncid,i,r8vals);call netcdf_err(error)
2236	mpp_file(unit)%time_values(1:len) = r8vals(1:len)
2237	case default
2238	call mpp_error( FATAL, 'Invalid data type for dimension' )
2239	end select
2240	endif
2241	! assign dimension atts
2242	if( nvatts.GT.0 )allocate(Axis(dimid)%Att(nvatts))
2243
2244	do j=1,nvatts
2245	Axis(dimid)%Att(j) = default_att
2246	enddo
2247
2248	do j=1,nvatts
2249	error=NF_INQ_ATTNAME(ncid,i,j,attname);call netcdf_err(error)
2250	error=NF_INQ_ATT(ncid,i,trim(attname),type,len);call netcdf_err(error)
2251
2252	Axis(dimid)%Att(j)%name = trim(attname)
2253	Axis(dimid)%Att(j)%type = type
2254	Axis(dimid)%Att(j)%len = len
2255
2256	select case (type)
2257	case (NF_CHAR)
2258	if (len.gt.512) call mpp_error(FATAL,'DIM ATT too long')
2259	error=NF_GET_ATT_TEXT(ncid,i,trim(attname),Axis(dimid)%Att(j)%catt);call netcdf_err(error)
2260	if( verbose .and. pe == 0 ) &
2261	print *, 'AXIS ',trim(Axis(dimid)%name),' ATT ',trim(attname),' ',Axis(dimid)%Att(j)%catt(1:len)
2262	! store integers in float arrays
2263	! assume dimension data not packed
2264	case (NF_SHORT)
2265	allocate(Axis(dimid)%Att(j)%fatt(len))
2266	error=NF_GET_ATT_INT2(ncid,i,trim(attname),i2vals);call netcdf_err(error)
2267	Axis(dimid)%Att(j)%fatt(1:len)=i2vals(1:len)
2268	if( verbose .and. pe == 0 ) &
2269	print *, 'AXIS ',trim(Axis(dimid)%name),' ATT ',trim(attname),' ',Axis(dimid)%Att(j)%fatt
2270	case (NF_INT)
2271	allocate(Axis(dimid)%Att(j)%fatt(len))
2272	error=NF_GET_ATT_INT(ncid,i,trim(attname),ivals);call netcdf_err(error)
2273	Axis(dimid)%Att(j)%fatt(1:len)=ivals(1:len)
2274	if( verbose .and. pe == 0 ) &
2275	print *, 'AXIS ',trim(Axis(dimid)%name),' ATT ',trim(attname),' ',Axis(dimid)%Att(j)%fatt
2276	case (NF_FLOAT)
2277	allocate(Axis(dimid)%Att(j)%fatt(len))
2278	error=NF_GET_ATT_REAL(ncid,i,trim(attname),rvals);call netcdf_err(error)
2279	Axis(dimid)%Att(j)%fatt(1:len)=rvals(1:len)
2280	if( verbose .and. pe == 0 ) &
2281	print *, 'AXIS ',trim(Axis(dimid)%name),' ATT ',trim(attname),' ',Axis(dimid)%Att(j)%fatt
2282	case (NF_DOUBLE)
2283	allocate(Axis(dimid)%Att(j)%fatt(len))
2284	error=NF_GET_ATT_DOUBLE(ncid,i,trim(attname),r8vals);call netcdf_err(error)
2285	Axis(dimid)%Att(j)%fatt(1:len)=r8vals(1:len)
2286	if( verbose .and. pe == 0 ) &
2287	print *, 'AXIS ',trim(Axis(dimid)%name),' ATT ',trim(attname),' ',Axis(dimid)%Att(j)%fatt
2288	case default
2289	call mpp_error( FATAL, 'Invalid data type for dimension at' )
2290	end select
2291	! assign pre-defined axis attributes
2292	select case(trim(attname))
2293	case('long_name')
2294	Axis(dimid)%longname=Axis(dimid)%Att(j)%catt(1:len)
2295	case('units')
2296	Axis(dimid)%units=Axis(dimid)%Att(j)%catt(1:len)
2297	case('cartesian_axis')
2298	Axis(dimid)%cartesian=Axis(dimid)%Att(j)%catt(1:len)
2299	case('positive')
2300	attval = Axis(dimid)%Att(j)%catt(1:len)
2301	if( attval.eq.'down' )then
2302	Axis(dimid)%sense=-1
2303	else if( attval.eq.'up' )then
2304	Axis(dimid)%sense=1
2305	endif
2306	end select
2307
2308	enddo
2309	! store axis info in filetype
2310	mpp_file(unit)%Axis(dimid) = Axis(dimid)
2311	endif
2312	enddo
2313	! assign variable info
2314	nv = 0
2315	do i=1, nvar_total
2316	error=NF_INQ_VAR(ncid,i,name,type,nvdims,dimids,nvatts);call netcdf_err(error)
2317	!
2318	! is this a dimension variable?
2319	!
2320	isdim=.false.
2321	do j=1,ndim
2322	if( trim(lowercase(name)).EQ.trim(lowercase(Axis(j)%name)) )isdim=.true.
2323	enddo
2324
2325	if( .not.isdim )then
2326	! for non-dimension variables
2327	nv=nv+1; if( nv.GT.mpp_file(unit)%nvar )call mpp_error( FATAL, 'variable index exceeds number of defined variables' )
2328	mpp_file(unit)%Var(nv)%type = type
2329	mpp_file(unit)%Var(nv)%id = i
2330	mpp_file(unit)%Var(nv)%name = name
2331	mpp_file(unit)%Var(nv)%natt = nvatts
2332	! determine packing attribute based on NetCDF variable type
2333	select case (type)
2334	case(NF_SHORT)
2335	mpp_file(unit)%Var(nv)%pack = 4
2336	case(NF_FLOAT)
2337	mpp_file(unit)%Var(nv)%pack = 2
2338	case(NF_DOUBLE)
2339	mpp_file(unit)%Var(nv)%pack = 1
2340	case (NF_INT)
2341	mpp_file(unit)%Var(nv)%pack = 2
2342	case default
2343	call mpp_error( FATAL, 'Invalid variable type in NetCDF file' )
2344	end select
2345	! assign dimension ids
2346	mpp_file(unit)%Var(nv)%ndim = nvdims
2347	allocate(mpp_file(unit)%Var(nv)%axes(nvdims))
2348	do j=1,nvdims
2349	mpp_file(unit)%Var(nv)%axes(j) = Axis(dimids(j))
2350	enddo
2351	allocate(mpp_file(unit)%Var(nv)%size(nvdims))
2352
2353	do j=1,nvdims
2354	if( dimids(j).eq.mpp_file(unit)%recdimid )then
2355	mpp_file(unit)%Var(nv)%time_axis_index = dimids(j)
2356	mpp_file(unit)%Var(nv)%size(j)=1 ! dimid length set to 1 here for consistency w/ mpp_write
2357	else
2358	mpp_file(unit)%Var(nv)%size(j)=Axis(dimids(j))%len
2359	endif
2360	enddo
2361	! assign variable atts
2362	if( nvatts.GT.0 )allocate(mpp_file(unit)%Var(nv)%Att(nvatts))
2363
2364	do j=1,nvatts
2365	mpp_file(unit)%Var(nv)%Att(j) = default_att
2366	enddo
2367
2368	do j=1,nvatts
2369	error=NF_INQ_ATTNAME(ncid,i,j,attname);call netcdf_err(error)
2370	error=NF_INQ_ATT(ncid,i,attname,type,len);call netcdf_err(error)
2371	mpp_file(unit)%Var(nv)%Att(j)%name = trim(attname)
2372	mpp_file(unit)%Var(nv)%Att(j)%type = type
2373	mpp_file(unit)%Var(nv)%Att(j)%len = len
2374
2375	select case (type)
2376	case (NF_CHAR)
2377	if (len.gt.512) call mpp_error(FATAL,'VAR ATT too long')
2378	error=NF_GET_ATT_TEXT(ncid,i,trim(attname),mpp_file(unit)%Var(nv)%Att(j)%catt(1:len));call netcdf_err(error)
2379	if (verbose .and. pe == 0 )&
2380	print *, 'Var ',nv,' ATT ',trim(attname),' ',mpp_file(unit)%Var(nv)%Att(j)%catt(1:len)
2381	! store integers as float internally
2382	case (NF_SHORT)
2383	allocate(mpp_file(unit)%Var(nv)%Att(j)%fatt(len))
2384	error=NF_GET_ATT_INT2(ncid,i,trim(attname),i2vals);call netcdf_err(error)
2385	mpp_file(unit)%Var(nv)%Att(j)%fatt(1:len)= i2vals(1:len)
2386	if( verbose .and. pe == 0 )&
2387	print *, 'Var ',nv,' ATT ',trim(attname),' ',mpp_file(unit)%Var(nv)%Att(j)%fatt
2388	case (NF_INT)
2389	allocate(mpp_file(unit)%Var(nv)%Att(j)%fatt(len))
2390	error=NF_GET_ATT_INT(ncid,i,trim(attname),ivals);call netcdf_err(error)
2391	mpp_file(unit)%Var(nv)%Att(j)%fatt(1:len)=ivals(1:len)
2392	if( verbose .and. pe == 0 )&
2393	print *, 'Var ',nv,' ATT ',trim(attname),' ',mpp_file(unit)%Var(nv)%Att(j)%fatt
2394	case (NF_FLOAT)
2395	allocate(mpp_file(unit)%Var(nv)%Att(j)%fatt(len))
2396	error=NF_GET_ATT_REAL(ncid,i,trim(attname),rvals);call netcdf_err(error)
2397	mpp_file(unit)%Var(nv)%Att(j)%fatt(1:len)=rvals(1:len)
2398	if( verbose .and. pe == 0 )&
2399	print *, 'Var ',nv,' ATT ',trim(attname),' ',mpp_file(unit)%Var(nv)%Att(j)%fatt
2400	case (NF_DOUBLE)
2401	allocate(mpp_file(unit)%Var(nv)%Att(j)%fatt(len))
2402	error=NF_GET_ATT_DOUBLE(ncid,i,trim(attname),r8vals);call netcdf_err(error)
2403	mpp_file(unit)%Var(nv)%Att(j)%fatt(1:len)=r8vals(1:len)
2404	if( verbose .and. pe == 0 ) &
2405	print *, 'Var ',nv,' ATT ',trim(attname),' ',mpp_file(unit)%Var(nv)%Att(j)%fatt
2406	case default
2407	call mpp_error( FATAL, 'Invalid data type for variable att' )
2408	end select
2409	! assign pre-defined field attributes
2410	select case (trim(attname))
2411	case ('long_name')
2412	mpp_file(unit)%Var(nv)%longname=mpp_file(unit)%Var(nv)%Att(j)%catt(1:len)
2413	case('units')
2414	mpp_file(unit)%Var(nv)%units=mpp_file(unit)%Var(nv)%Att(j)%catt(1:len)
2415	case('scale_factor')
2416	mpp_file(unit)%Var(nv)%scale=mpp_file(unit)%Var(nv)%Att(j)%fatt(1)
2417	case('missing')
2418	mpp_file(unit)%Var(nv)%missing=mpp_file(unit)%Var(nv)%Att(j)%fatt(1)
2419	case('add_offset')
2420	mpp_file(unit)%Var(nv)%add=mpp_file(unit)%Var(nv)%Att(j)%fatt(1)
2421	case('valid_range')
2422	mpp_file(unit)%Var(nv)%min=mpp_file(unit)%Var(nv)%Att(j)%fatt(1)
2423	mpp_file(unit)%Var(nv)%max=mpp_file(unit)%Var(nv)%Att(j)%fatt(2)
2424	end select
2425	enddo
2426	endif
2427	enddo ! end variable loop
2428	else
2429	call mpp_error( FATAL, 'MPP READ CURRENTLY DOES NOT SUPPORT NON-NETCDF' )
2430	endif
2431
2432	mpp_file(unit)%initialized = .TRUE.
2433	#else
2434	call mpp_error( FATAL, 'MPP_READ currently requires use_netCDF option' )
2435	#endif
2436	return
2437	end subroutine mpp_read_meta
2438
2439
2440	subroutine mpp_get_info( unit, ndim, nvar, natt, ntime )
2441
2442	integer, intent(in) :: unit
2443	integer, intent(out) :: ndim, nvar, natt, ntime
2444
2445
2446	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_GET_INFO: must first call mpp_io_init.' )
2447	if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_GET_INFO: invalid unit number.' )
2448
2449	ndim = mpp_file(unit)%ndim
2450	nvar = mpp_file(unit)%nvar
2451	natt = mpp_file(unit)%natt
2452	ntime = mpp_file(unit)%time_level
2453
2454	return
2455
2456	end subroutine mpp_get_info
2457
2458
2459	subroutine mpp_get_global_atts( unit, global_atts )
2460	!
2461	! copy global file attributes for use by user
2462	!
2463	! global_atts is an attribute type which is allocated from the
2464	! calling routine
2465
2466	integer, intent(in) :: unit
2467	type(atttype), intent(inout) :: global_atts(:)
2468	integer :: natt,i
2469
2470	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_GET_INFO: must first call mpp_io_init.' )
2471	if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_GET_INFO: invalid unit number.' )
2472
2473	if (size(global_atts).lt.mpp_file(unit)%natt) &
2474	call mpp_error(FATAL, 'MPP_GET_ATTS: atttype not dimensioned properly in calling routine')
2475
2476	natt = mpp_file(unit)%natt
2477	global_atts = default_att
2478
2479	do i=1,natt
2480	global_atts(i) = mpp_file(unit)%Att(i)
2481	enddo
2482
2483	return
2484	end subroutine mpp_get_global_atts
2485
2486	subroutine mpp_get_field_atts( field, name, units, longname, min, max, missing, ndim, siz, axes, atts )
2487
2488	type(fieldtype), intent(in) :: field
2489	character(len=*), intent(out) , optional :: name, units
2490	character(len=*), intent(out), optional :: longname
2491	real,intent(out), optional :: min,max,missing
2492	integer, intent(out), optional :: ndim
2493	integer, intent(out), dimension(:), optional :: siz
2494
2495	type(atttype), intent(out), optional, dimension(:) :: atts
2496	type(axistype), intent(out), optional, dimension(:) :: axes
2497
2498	integer :: n,m
2499
2500	if (PRESENT(name)) name = field%name
2501	if (PRESENT(units)) units = field%units
2502	if (PRESENT(longname)) longname = field%longname
2503	if (PRESENT(min)) min = field%min
2504	if (PRESENT(max)) max = field%max
2505	if (PRESENT(missing)) missing = field%missing
2506	if (PRESENT(ndim)) ndim = field%ndim
2507	if (PRESENT(atts)) then
2508	atts = default_att
2509	n = size(atts);m=size(field%Att)
2510	if (n.LT.m) call mpp_error(FATAL,'attribute array not large enough in mpp_get_field_atts')
2511	atts(1:m) = field%Att(1:m)
2512	end if
2513	if (PRESENT(axes)) then
2514	axes = default_axis
2515	n = size(axes);m=field%ndim
2516	if (n.LT.m) call mpp_error(FATAL,'axis array not large enough in mpp_get_field_atts')
2517	axes(1:m) = field%axes(1:m)
2518	end if
2519	if (PRESENT(siz)) then
2520	siz = -1
2521	n = size(siz);m=field%ndim
2522	if (n.LT.m) call mpp_error(FATAL,'size array not large enough in mpp_get_field_atts')
2523	siz(1:m) = field%size(1:m)
2524	end if
2525	return
2526	end subroutine mpp_get_field_atts
2527
2528	subroutine mpp_get_axis_atts( axis, name, units, longname, cartesian, sense, len, natts, atts )
2529
2530	type(axistype), intent(in) :: axis
2531	character(len=*), intent(out) , optional :: name, units
2532	character(len=*), intent(out), optional :: longname, cartesian
2533	integer,intent(out), optional :: sense, len , natts
2534	type(atttype), intent(out), optional, dimension(:) :: atts
2535
2536	integer :: n,m
2537
2538	if (PRESENT(name)) name = axis%name
2539	if (PRESENT(units)) units = axis%units
2540	if (PRESENT(longname)) longname = axis%longname
2541	if (PRESENT(cartesian)) cartesian = axis%cartesian
2542	if (PRESENT(sense)) sense = axis%sense
2543	if (PRESENT(len)) len = axis%len
2544	if (PRESENT(atts)) then
2545	atts = default_att
2546	n = size(atts);m=size(axis%Att)
2547	if (n.LT.m) call mpp_error(FATAL,'attribute array not large enough in mpp_get_field_atts')
2548	atts(1:m) = axis%Att(1:m)
2549	end if
2550	if (PRESENT(natts)) natts = size(axis%Att)
2551
2552	return
2553	end subroutine mpp_get_axis_atts
2554
2555
2556	subroutine mpp_get_fields( unit, variables )
2557	!
2558	! copy variable information from file (excluding data)
2559	! global_atts is an attribute type which is allocated from the
2560	! calling routine
2561	!
2562	integer, intent(in) :: unit
2563	type(fieldtype), intent(inout) :: variables(:)
2564
2565	integer :: nvar,i
2566
2567	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_GET_FIELDS: must first call mpp_io_init.' )
2568	if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_GET_FIELDS: invalid unit number.' )
2569
2570	if (size(variables).ne.mpp_file(unit)%nvar) &
2571	call mpp_error(FATAL, 'MPP_GET_FIELDS: fieldtype not dimensioned properly in calling routine')
2572
2573	nvar = mpp_file(unit)%nvar
2574
2575	do i=1,nvar
2576	variables(i) = mpp_file(unit)%Var(i)
2577	enddo
2578
2579	return
2580	end subroutine mpp_get_fields
2581
2582	subroutine mpp_get_axes( unit, axes, time_axis )
2583	!
2584	! copy variable information from file (excluding data)
2585	! global_atts is an attribute type which is allocated from the
2586	! calling routine
2587	!
2588	integer, intent(in) :: unit
2589	type(axistype), intent(out) :: axes(:)
2590	type(axistype), intent(out), optional :: time_axis
2591	character(len=128) :: name
2592	logical :: save
2593	integer :: ndim,i, nvar, j, num_dims, k
2594
2595	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_GET_AXES: must first call mpp_io_init.' )
2596	if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_GET_AXES: invalid unit number.' )
2597
2598	if (size(axes).ne.mpp_file(unit)%ndim) &
2599	call mpp_error(FATAL, 'MPP_GET_AXES: axistype not dimensioned properly in calling routine')
2600
2601
2602	if (PRESENT(time_axis)) time_axis = default_axis
2603	ndim = mpp_file(unit)%ndim
2604	do i=1,ndim
2605	if (ASSOCIATED(mpp_file(unit)%Axis(i)%data)) then
2606	axes(i)=mpp_file(unit)%Axis(i)
2607	else
2608	axes(i)=mpp_file(unit)%Axis(i)
2609	if (PRESENT(time_axis)) time_axis = mpp_file(unit)%Axis(i)
2610	endif
2611	enddo
2612
2613	return
2614	end subroutine mpp_get_axes
2615
2616	subroutine mpp_get_times( unit, time_values )
2617	!
2618	! copy time information from file and convert to time_type
2619	!
2620	integer, intent(in) :: unit
2621	real(DOUBLE_KIND), intent(inout) :: time_values(:)
2622
2623	integer :: ntime,i
2624
2625	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_GET_TIMES: must first call mpp_io_init.' )
2626	if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_GET_TIMES: invalid unit number.' )
2627
2628	if (size(time_values).ne.mpp_file(unit)%time_level) &
2629	call mpp_error(FATAL, 'MPP_GET_TIMES: time_values not dimensioned properly in calling routine')
2630
2631	ntime = mpp_file(unit)%time_level
2632
2633	do i=1,ntime
2634	time_values(i) = mpp_file(unit)%time_values(i)
2635	enddo
2636
2637
2638
2639	return
2640	end subroutine mpp_get_times
2641
2642	function mpp_get_field_index(fields,fieldname)
2643
2644	type(fieldtype), dimension(:) :: fields
2645	character(len=*) :: fieldname
2646	integer :: mpp_get_field_index
2647
2648	integer :: n
2649
2650	mpp_get_field_index = -1
2651
2652	do n=1,size(fields)
2653	if (lowercase(fields(n)%name) == lowercase(fieldname)) then
2654	mpp_get_field_index = n
2655	exit
2656	endif
2657	enddo
2658
2659	return
2660	end function mpp_get_field_index
2661
2662	function mpp_get_field_size(field)
2663
2664	type(fieldtype) :: field
2665	integer :: mpp_get_field_size(4)
2666
2667	integer :: n
2668
2669	mpp_get_field_size = -1
2670
2671	mpp_get_field_size(1) = field%size(1)
2672	mpp_get_field_size(2) = field%size(2)
2673	mpp_get_field_size(3) = field%size(3)
2674	mpp_get_field_size(4) = field%size(4)
2675
2676	return
2677	end function mpp_get_field_size
2678
2679
2680	subroutine mpp_get_axis_data( axis, data )
2681
2682	type(axistype), intent(in) :: axis
2683	real, dimension(:), intent(out) :: data
2684
2685
2686	if (size(data).lt.axis%len) call mpp_error(FATAL,'MPP_GET_AXIS_DATA: data array not large enough')
2687	if (.NOT.ASSOCIATED(axis%data)) then
2688	call mpp_error(NOTE,'MPP_GET_AXIS_DATA: use mpp_get_times for record dims')
2689	data = 0.
2690	else
2691	data(1:axis%len) = axis%data
2692	endif
2693
2694	return
2695	end subroutine mpp_get_axis_data
2696
2697
2698	function mpp_get_recdimid(unit)
2699	!
2700	integer, intent(in) :: unit
2701	integer :: mpp_get_recdimid
2702
2703
2704	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_GET_RECDIMID: must first call mpp_io_init.' )
2705	if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_GET_RECDIMID: invalid unit number.' )
2706
2707	mpp_get_recdimid = mpp_file(unit)%recdimid
2708
2709	return
2710	end function mpp_get_recdimid
2711
2712	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2713	! !
2714	! mpp_get_iospec, mpp_flush: OS-dependent calls !
2715	! !
2716	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2717
2718	subroutine mpp_flush(unit)
2719	!flush the output on a unit, syncing with disk
2720	integer, intent(in) :: unit
2721
2722	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_FLUSH: must first call mpp_io_init.' )
2723	if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_FLUSH: invalid unit number.' )
2724	if( .NOT.mpp_file(unit)%initialized )call mpp_error( FATAL, 'MPP_FLUSH: cannot flush a file during writing of metadata.' )
2725	if( mpp_file(unit)%threading.EQ.MPP_SINGLE .AND. pe.NE.mpp_root_pe() )return
2726
2727	if( mpp_file(unit)%format.EQ.MPP_NETCDF )then
2728	#ifdef use_netCDF
2729	error = NF_SYNC(mpp_file(unit)%ncid); call netcdf_err(error)
2730	#endif
2731	else
2732	call FLUSH(unit)
2733	end if
2734	return
2735	end subroutine mpp_flush
2736
2737	subroutine mpp_get_iospec( unit, iospec )
2738	integer, intent(in) :: unit
2739	character(len=*), intent(out) :: iospec
2740
2741	if( .NOT.module_is_initialized )call mpp_error( FATAL, 'MPP_GET_IOSPEC: must first call mpp_io_init.' )
2742	if( .NOT.mpp_file(unit)%opened )call mpp_error( FATAL, 'MPP_GET_IOSPEC: invalid unit number.' )
2743	#ifdef SGICRAY
2744	!currently will write to stdout: don't know how to trap and return as string to iospec
2745	call ASSIGN( 'assign -V f:'//trim(mpp_file(unit)%name), error )
2746	#endif
2747	return
2748	end subroutine mpp_get_iospec
2749
2750	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2751	! !
2752	! netCDF-specific routines: mpp_get_id, netcdf_error !
2753	! !
2754	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2755
2756	function mpp_get_ncid(unit)
2757	integer :: mpp_get_ncid
2758	integer, intent(in) :: unit
2759
2760	mpp_get_ncid = mpp_file(unit)%ncid
2761	return
2762	end function mpp_get_ncid
2763
2764	function mpp_get_axis_id(axis)
2765	integer mpp_get_axis_id
2766	type(axistype), intent(in) :: axis
2767	mpp_get_axis_id = axis%id
2768	return
2769	end function mpp_get_axis_id
2770
2771	function mpp_get_field_id(field)
2772	integer mpp_get_field_id
2773	type(fieldtype), intent(in) :: field
2774	mpp_get_field_id = field%id
2775	return
2776	end function mpp_get_field_id
2777
2778	subroutine netcdf_err(err)
2779	integer, intent(in) :: err
2780	character(len=80) :: errmsg
2781	integer :: unit
2782
2783	#ifdef use_netCDF
2784	if( err.EQ.NF_NOERR )return
2785	errmsg = NF_STRERROR(err)
2786	call mpp_io_exit() !make sure you close all open files
2787	call mpp_error( FATAL, 'NETCDF ERROR: '//trim(errmsg) )
2788	#endif
2789	return
2790	end subroutine netcdf_err
2791
2792	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2793	! !
2794	! minor routines: mpp_get_unit_range, mpp_set_unit_range !
2795	! !
2796	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2797
2798	subroutine mpp_get_unit_range( unit_begin_out, unit_end_out )
2799	integer, intent(out) :: unit_begin_out, unit_end_out
2800
2801	unit_begin_out = unit_begin; unit_end_out = unit_end
2802	return
2803	end subroutine mpp_get_unit_range
2804
2805	subroutine mpp_set_unit_range( unit_begin_in, unit_end_in )
2806	integer, intent(in) :: unit_begin_in, unit_end_in
2807
2808	if( unit_begin_in.GT.unit_end_in )call mpp_error( FATAL, 'MPP_SET_UNIT_RANGE: unit_begin_in.GT.unit_end_in.' )
2809	if( unit_begin_in.LT.0 )call mpp_error( FATAL, 'MPP_SET_UNIT_RANGE: unit_begin_in.LT.0.' )
2810	if( unit_end_in .GT.maxunits )call mpp_error( FATAL, 'MPP_SET_UNIT_RANGE: unit_end_in.GT.maxunits.' )
2811	unit_begin = unit_begin_in; unit_end = unit_end_in
2812	return
2813	end subroutine mpp_set_unit_range
2814
2815	subroutine mpp_modify_axis_meta( axis, name, units, longname, cartesian, data )
2816
2817	type(axistype), intent(inout) :: axis
2818	character(len=*), intent(in), optional :: name, units, longname, cartesian
2819	real, dimension(:), intent(in), optional :: data
2820
2821	if (PRESENT(name)) axis%name = trim(name)
2822	if (PRESENT(units)) axis%units = trim(units)
2823	if (PRESENT(longname)) axis%longname = trim(longname)
2824	if (PRESENT(cartesian)) axis%cartesian = trim(cartesian)
2825	if (PRESENT(data)) then
2826	axis%len = size(data)
2827	if (ASSOCIATED(axis%data)) deallocate(axis%data)
2828	allocate(axis%data(axis%len))
2829	axis%data = data
2830	endif
2831
2832	return
2833	end subroutine mpp_modify_axis_meta
2834
2835	subroutine mpp_modify_field_meta( field, name, units, longname, min, max, missing, axes )
2836
2837	type(fieldtype), intent(inout) :: field
2838	character(len=*), intent(in), optional :: name, units, longname
2839	real, intent(in), optional :: min, max, missing
2840	type(axistype), dimension(:), intent(inout), optional :: axes
2841
2842	if (PRESENT(name)) field%name = trim(name)
2843	if (PRESENT(units)) field%units = trim(units)
2844	if (PRESENT(longname)) field%longname = trim(longname)
2845	if (PRESENT(min)) field%min = min
2846	if (PRESENT(max)) field%max = max
2847	if (PRESENT(missing)) field%missing = missing
2848	! if (PRESENT(axes)) then
2849	! axis%len = size(data)
2850	! deallocate(axis%data)
2851	! allocate(axis%data(axis%len))
2852	! axis%data = data
2853	! endif
2854
2855	return
2856	end subroutine mpp_modify_field_meta
2857
2858	function lowercase (cs)
2859	character(len=*), intent(in) :: cs
2860	character(len=len(cs)) :: lowercase
2861	character :: ca(len(cs))
2862
2863	integer, parameter :: co=iachar('a')-iachar('A') ! case offset
2864
2865	ca = transfer(cs,"x",len(cs))
2866	where (ca >= "A" .and. ca <= "Z") ca = achar(iachar(ca)+co)
2867	lowercase = transfer(ca,cs)
2868
2869	end function lowercase
2870
2871	end module mpp_io_mod
2872
2873	#ifdef test_mpp_io
2874	program mpp_io_test
2875
2876	use mpp_mod
2877	use mpp_domains_mod
2878	use mpp_io_mod
2879
2880	implicit none
2881	#include <mpif.h>
2882
2883	#ifdef use_netCDF
2884	#include <netcdf.inc>
2885	#endif
2886
2887	integer :: pe, npes
2888	type(domain2D) :: domain
2889	integer :: nx=240, ny=240, nz=40, nt=2, halo=2, stackmax=32768, stackmaxd=32768
2890	real, dimension(:,:,:), allocatable :: data, gdata, rdata
2891	real, dimension(:,:), allocatable :: lat,lon
2892	integer :: is, ie, js, je, isd, ied, jsd, jed
2893	integer :: tk, tk0, tks_per_sec
2894	integer :: i,j,k, unit=7, layout(2)
2895	logical :: debug=.FALSE., opened
2896	character(len=64) :: file='test', iospec='-F cachea', varname
2897	namelist / mpp_io_nml / nx, ny, nz, nt, halo, stackmax, stackmaxd, debug, file, iospec
2898	integer :: ndim, nvar, natt, ntime
2899	type(atttype), allocatable :: atts(:)
2900	type(fieldtype), allocatable :: vars(:)
2901	type(axistype), allocatable :: axes(:)
2902	real(DOUBLE_KIND), allocatable :: tstamp(:)
2903	real(DOUBLE_KIND) :: time
2904	real :: dt,bytes,bw
2905	type(axistype) :: x, y, z, t
2906	type(fieldtype) :: f,latij,lonij
2907	type(domain1D) :: xdom, ydom
2908	integer(LONG_KIND) :: rchk, chk
2909
2910	call mpi_init(is)
2911	call mpp_init(mpp_comm=MPI_COMM_WORLD)
2912	pe = mpp_pe()
2913	npes = mpp_npes()
2914
2915	!possibly open a file called mpp_io.nml
2916	do
2917	inquire( unit=unit, opened=opened )
2918	if( .NOT.opened )exit
2919	unit = unit + 1
2920	if( unit.EQ.100 )call mpp_error( FATAL, 'Unable to locate unit number.' )
2921	end do
2922	open( unit=unit, status='OLD', file='mpp_io.nml', err=10 )
2923	read( unit,mpp_io_nml )
2924	close(unit)
2925	10 continue
2926
2927	call SYSTEM_CLOCK( count_rate=tks_per_sec )
2928	if( debug )then
2929	call mpp_io_init(MPP_DEBUG)
2930	else
2931	call mpp_io_init()
2932	end if
2933	call mpp_set_stack_size(stackmax)
2934	call mpp_domains_set_stack_size(stackmaxd)
2935
2936	if( pe.EQ.mpp_root_pe() )then
2937	print '(a,6i4)', 'npes, nx, ny, nz, nt, halo=', npes, nx, ny, nz, nt, halo
2938	print *, 'Using NEW domaintypes and calls...'
2939	end if
2940	!define global data array
2941	allocate( gdata(nx,ny,nz) )
2942	if( pe.EQ.mpp_root_pe() )then
2943	! call random_number(gdata) )
2944	!fill in global array: with k.iiijjj
2945	gdata = 0.
2946	do k = 1,nz
2947	do j = 1,ny
2948	do i = 1,nx
2949	gdata(i,j,k) = k + i1e-3 + j1e-6
2950	end do
2951	end do
2952	end do
2953	end if
2954	call mpp_broadcast( gdata, size(gdata), mpp_root_pe() )
2955
2956	!define domain decomposition
2957	call mpp_define_layout( (/1,nx,1,ny/), npes, layout )
2958	call mpp_define_domains( (/1,nx,1,ny/), layout, domain, xhalo=halo, yhalo=halo )
2959	call mpp_get_compute_domain( domain, is, ie, js, je )
2960	call mpp_get_data_domain ( domain, isd, ied, jsd, jed )
2961	call mpp_get_domain_components( domain, xdom, ydom )
2962	allocate( data(isd:ied,jsd:jed,nz) )
2963	data(is:ie,js:je,:) = gdata(is:ie,js:je,:)
2964
2965	allocate( lat(isd:ied,jsd:jed) )
2966	allocate( lon(isd:ied,jsd:jed) )
2967	do j=js,je
2968	do i=is,ie
2969	lat(i,j)=i+j
2970	lon(i,j)=-(i+j)
2971	enddo
2972	enddo
2973
2974	!tests
2975	write( file,'(a,i3.3)' )trim(file), npes
2976
2977
2978	!netCDF distributed write
2979	if( pe.EQ.mpp_root_pe() )print *, 'netCDF distributed write'
2980	call mpp_open( unit, trim(file)//'d', action=MPP_OVERWR, form=MPP_NETCDF, threading=MPP_MULTI, fileset=MPP_MULTI )
2981	call mpp_write_meta( unit, x,'X','km','x-coordinate in Cartesian system',&
2982	domain=xdom, data=(/(i-1.,i=1,nx)/) )
2983	call mpp_write_meta( unit, y,'Y','km','x-coordinate in Cartesian system',&
2984	domain=ydom, data=(/(i-1.,i=1,ny)/) )
2985
2986	call mpp_write_meta( unit, z,'Z','km','height in Cartesian system ',&
2987	sense=1,data=(/(i-1.,i=1,nz)/) )
2988	! The optional parameter sense=1 has the same meaning of writing
2989	! call mpp_write_meta( unit, z%id, 'positive', cval='up')
2990
2991	call mpp_write_meta( unit, t, 'time', 'sec', 'Time' )
2992
2993	call mpp_write_meta( unit, latij, (/x,y/), 'lat', 'degrees_north', &
2994	'latitude')
2995	call mpp_write_meta( unit, lonij, (/x,y/), 'lon', 'degrees_east', &
2996	'longitude')
2997
2998	call mpp_write_meta( unit, f, (/x,y,z,t/), 'T', 'K', 'Temperature', pack=1 )
2999	call mpp_write_meta( unit, mpp_get_id(f), 'coordinates',cval='lon lat')
3000	call mpp_write( unit, x )
3001	call mpp_write( unit, y )
3002	call mpp_write( unit, z )
3003	call mpp_write( unit, latij, domain, lat)
3004	call mpp_write( unit, lonij, domain, lon)
3005	call system_clock(tk0)
3006	do i = 0,nt-1
3007	time = i*10.
3008	call mpp_write( unit, f, domain, data, time )
3009	end do
3010	call mpp_close(unit)
3011	call system_clock(tk)
3012	dt = real(tk-tk0)/(tks_per_sec)
3013	dt = max( dt, epsilon(dt) )
3014	bytes=(ie-is+1)(je-js+1)nznpesnt*kind(data)
3015	bw=bytes*1.e-6/dt
3016	if( pe.EQ.mpp_root_pe() ) &
3017	write(stdout(),'(a,g12.5,a,g12.5,a,g12.5,a)')'MPP_WRITE multi: n = ',&
3018	bytes,' bytes time =',dt,' s bw = ',bw,' MByte/s'
3019
3020	!netCDF single-threaded write
3021	if( pe.EQ.mpp_root_pe() )print *, 'netCDF single-threaded write'
3022	call mpp_open( unit, trim(file)//'s', action=MPP_OVERWR, form=MPP_NETCDF, threading=MPP_SINGLE )
3023	! call mpp_write_meta( unit, x, 'X', 'km', 'X distance', domain=xdom, data=(/(i-1.,i=1,nx)/) )
3024	! call mpp_write_meta( unit, y, 'Y', 'km', 'Y distance', domain=ydom, data=(/(i-1.,i=1,ny)/) )
3025	! call mpp_write_meta( unit, z, 'Z', 'km', 'Z distance', data=(/(i-1.,i=1,nz)/) )
3026	! call mpp_write_meta( unit, t, 'T', 'sec', 'Time' )
3027	! call mpp_write_meta( unit, f, (/x,y,z,t/), 'Data', 'metres', 'Random data', pack=1 )
3028	! call mpp_write( unit, x )
3029	! call mpp_write( unit, y )
3030	! call mpp_write( unit, z )
3031	call mpp_write_meta( unit, x,'X','km','x-coordinate in Cartesian system',&
3032	domain=xdom, data=(/(i-1.,i=1,nx)/) )
3033	call mpp_write_meta( unit, y,'Y','km','x-coordinate in Cartesian system',&
3034	domain=ydom, data=(/(i-1.,i=1,ny)/) )
3035
3036	call mpp_write_meta( unit, z,'Z','km','height in Cartesian system ',&
3037	sense=1,data=(/(i-1.,i=1,nz)/) )
3038
3039	! The optional parameter sense=1 has the same meaning as writing
3040	! call mpp_write_meta( unit, z%id, 'positive', cval='up')
3041
3042	call mpp_write_meta( unit, t, 'time', 'sec', 'Time' )
3043
3044	call mpp_write_meta( unit, latij, (/x,y/), 'lat', 'degrees_north', &
3045	'latitude')
3046	call mpp_write_meta( unit, lonij, (/x,y/), 'lon', 'degrees_east', &
3047	'longitude')
3048
3049	call mpp_write_meta( unit, f, (/x,y,z,t/), 'T', 'K', 'Temperature', pack=1 )
3050	call mpp_write_meta( unit, mpp_get_id(f), 'coordinates',cval='lon lat')
3051	call mpp_write( unit, x )
3052	call mpp_write( unit, y )
3053	call mpp_write( unit, z )
3054	call mpp_write( unit, latij, domain, lat)
3055	call mpp_write( unit, lonij, domain, lon)
3056
3057	call system_clock(tk0)
3058	do i = 0,nt-1
3059	time = i*10.
3060	call mpp_write( unit, f, domain, data, time )
3061	end do
3062	call mpp_close(unit)
3063	call system_clock(tk)
3064	dt = real(tk-tk0)/(tks_per_sec)
3065	dt = max( dt, epsilon(dt) )
3066	bytes=(ie-is+1)(je-js+1)nznpesnt*kind(data)
3067	bw=bytes*1.e-6/dt
3068	if( pe.EQ.mpp_root_pe() ) &
3069	write(stdout(),'(a,g12.5,a,g12.5,a,g12.5,a)')'MPP_WRITE single: n = ',&
3070	bytes,' bytes time =',dt,' s bw = ',bw,' MByte/s'
3071
3072	!netCDF multi-threaded read
3073	if( pe.EQ.mpp_root_pe() )print *, 'netCDF multi-threaded read'
3074	call mpp_sync() !wait for previous write to complete
3075	call mpp_open( unit, trim(file)//'s', action=MPP_RDONLY, form=MPP_NETCDF, threading=MPP_MULTI, fileset=MPP_SINGLE )
3076	call mpp_get_info( unit, ndim, nvar, natt, ntime )
3077	allocate( atts(natt) )
3078	allocate( axes(ndim) )
3079	allocate( vars(nvar) )
3080	allocate( tstamp(ntime) )
3081	call mpp_get_atts ( unit, atts(:) )
3082	call mpp_get_axes ( unit, axes(:) )
3083	call mpp_get_fields ( unit, vars(:) )
3084	call mpp_get_times( unit, tstamp(:) )
3085
3086	k=-1
3087	do i=1,nvar
3088	call mpp_get_atts(vars(i),name=varname)
3089	write(stdout(),*) varname
3090	if(varname.eq.'T')k=i
3091	enddo
3092
3093	! if( varname.NE.'Data' )call mpp_error( FATAL, 'File being read is not the expected one.' )
3094	if( k.lt.0 )call mpp_error( FATAL, 'File being read is not the expected one.' )
3095	allocate( rdata(is:ie,js:je,nz) )
3096	call system_clock(tk0)
3097	call mpp_read( unit, vars(k), domain, rdata, 1 )
3098	call system_clock(tk)
3099	dt = real(tk-tk0)/(tks_per_sec)
3100	dt = max( dt, epsilon(dt) )
3101	bytes=(ie-is+1)(je-js+1)nznpeskind(data)
3102	bw=bytes*1.e-6/dt
3103	if( pe.EQ.mpp_root_pe() ) &
3104	write(stdout(),'(a,g12.5,a,g12.5,a,g12.5,a)')'MPP_READ multi: n = ',&
3105	bytes,' bytes time =',dt,' s bw = ',bw,' MByte/s'
3106	rchk = mpp_chksum(rdata(is:ie,js:je,:))
3107	chk = mpp_chksum( data(is:ie,js:je,:))
3108	if( pe.EQ.mpp_root_pe() )print '(a,2z18)', 'checksum=', rchk, chk
3109	if( rchk.NE.chk )call mpp_error( FATAL, 'Checksum error on multi-threaded netCDF read.' )
3110
3111	call mpp_io_exit()
3112	call mpp_domains_exit()
3113	call mpp_exit()
3114	call mpi_finalize(is)
3115
3116	end program mpp_io_test
3117
3118	#endif

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: CPL/oasis3/trunk/src/lib/mpp_io/src/mpp_io_mod.F90 @ 1677

Download in other formats: